[MOAB-dev] AdaptiveKDTree and point search

[Vijay S. Mahadevan]

Lukasz, Thanks for verifying this on both versions (5.0, 5.0.1). For
use in DMMoab, can modify MOAB.py with the appropriate hash for
v5.0.1, or use your custom compiled/installed MOAB version, if you
don't use the --download-moab method. I am not sure whether its worth
it to apply a backfix to 5.0, but we could release a 5.0.0-p1 if
necessary, which still would need a hash update in PETSc/MOAB.py
anyway.

Vijay

On Wed, Nov 15, 2017 at 12:54 PM, Lukasz Kaczmarczyk
<Lukasz.Kaczmarczyk at glasgow.ac.uk> wrote:
> Vijay,
>
> Ok, I found the problem, which is fixed in MOAB 5.0.1 but not in MOAB 5.0.0 which goes with petsc installation. This is code snippet which you can use to reproduce the error. Problem is with HO nodes; the evaluator should ask for topological nodes only. Anyway, it works for linear geometry only. Let me know if you need mesh, but it will go with any tet mesh.
>
> No problem at the end, since is fine with MOAB 5.0.1.
>
> #include <moab/ElemEvaluator.hpp>
> #include <moab/SpatialLocator.hpp>
>
> int main(int argc, char *argv[]) {
>
>   moab::Core mb_instance;
>   moab::Interface &moab = mb_instance;
>
>   const char *option;
>   option = "";
>
>   moab.load_file("examples/sslv116/init_stress_tet/out.h5m", 0, option);
>
>   Range tets;
>   moab.get_entities_by_dimension(0, 3, tets);
>
>   EntityHandle meshset;
>   moab.create_meshset(MESHSET_SET, meshset);
>   moab.add_entities(meshset, tets);
>   moab.convert_entities(meshset, false, false, true);
>   moab.delete_entities(&meshset, 1);
>
>   AdaptiveKDTree tree(&moab);
>   cout << "Build tree ... ";
>   EntityHandle root;
>   tree.build_tree(tets, &root);
>   ElemEvaluator elem_eval(&moab);
>   tree.set_eval(&elem_eval);
>   cout << "done" << endl;
>
>   double coords[] = {-10, 0, 10};
>   CartVect params(coords);
>   EntityHandle leaf_out;
>   tree.point_search(coords, leaf_out, 3.87121 * 1e-5, 3.87121*1e-3, NULL, NULL,
>                     &params);
>
>   if (moab.type_from_handle(leaf_out) == MBTET) {
>     cout << "Tet found" << endl;
>   }
>
>   return 0;
> }
>
> Regards,
> Lukasz
>
>
>
>
>> On 15 Nov 2017, at 15:55, Vijay S. Mahadevan <vijay.m at gmail.com> wrote:
>>
>> Lukasz,
>>
>> Any chance you can isolate it to a small test ? Looks like the LD code
>> for LinearTet is hitting an assertion. We need to see why would nverts
>> for a LinearTet element not be 4. What kind of mesh are you using ?
>>
>> We've been refining and adding many new capabilities in the coupler
>> and LD recently. I just want to make sure that the issue below is not
>> something in the setup of the mesh, but a regression from the new
>> changes before we spend a lot of time finding the exact problem here.
>>
>> Vijay
>>
>> On Wed, Nov 15, 2017 at 10:36 AM, Lukasz Kaczmarczyk
>> <Lukasz.Kaczmarczyk at glasgow.ac.uk> wrote:
>>> Hello,
>>>
>>> I switch recently to MOAB 5.0.0 and I encounter some problems with point search.  This is for tet mesh, I expecting that I doing something wrong. It was working with older MOAB. I can work with leaf, but will be nice to have this ok.
>>>
>>>   myTree = boost::shared_ptr<AdaptiveKDTree>(new AdaptiveKDTree(&mwlsMoab));
>>>   PetscPrintf(mField.get_comm(),"Build tree ... ");
>>>   rval = myTree->build_tree(levelTets,&treeRoot); CHKERRQ_MOAB(rval);
>>>
>>>   boost::shared_ptr<ElemEvaluator> elemEval(myTree,new ElemEvaluator(&mwlsMoab));
>>>   myTree->set_eval(elemEval.get());
>>>
>>>   PetscPrintf(mField.get_comm(),"done\n”);
>>>
>>> and using it like that
>>>
>>> rval = myTree->point_search(
>>>         material_coords,leaf_out,maxEdgeL*1e-5,maxEdgeL*1e-3,NULL,NULL,&params
>>>       ); CHKERRQ_MOAB(rval);
>>>
>>> I get following error,
>>> ssertion failed: (nverts == 4 && verts), function initFcn, file LocalDiscretization/LinearTet.cpp, line 21.
>>> [likask:32071] *** Process received signal ***
>>> [likask:32071] Signal: Abort trap: 6 (6)
>>> [likask:32071] Signal code:  (0)
>>> [likask:32071] [ 0] 0   libsystem_platform.dylib            0x00007fff9d14f52a _sigtramp + 26
>>> [likask:32071] [ 1] 0   ???                                 0x00007fff5acbdea0 0x0 + 140734716698272
>>> [likask:32071] [ 2] 0   libsystem_c.dylib                   0x00007fff8aca76df abort + 129
>>> [likask:32071] [ 3] 0   libsystem_c.dylib                   0x00007fff8ac6edd8 basename + 0
>>> [likask:32071] [ 4] 0   libMOAB.0.dylib                     0x000000010a4fc1ec _ZN4moab9LinearTet7initFcnEPKdiRPd + 108
>>> [likask:32071] [ 5] 0   libMOAB.0.dylib                     0x000000010a1f5cad _ZN4moab13ElemEvaluator14set_ent_handleEm + 1309
>>> [likask:32071] [ 6] 0   libMOAB.0.dylib                     0x000000010a4f6e45 _ZN4moab13ElemEvaluator22find_containing_entityERNS_5RangeEPKdddRmPdPj + 197
>>> [likask:32071] [ 7] 0   libMOAB.0.dylib                     0x000000010a19c033 _ZN4moab13ElemEvaluator22find_containing_entityEmPKdddRmPdPj + 371
>>> [likask:32071] [ 8] 0   libMOAB.0.dylib                     0x000000010a18ef65 _ZN4moab14AdaptiveKDTree12point_searchEPKdRmddPbPmPNS_8CartVectE + 2069
>>> [likask:32071] [ 9] 0   mwls                                0x00000001051a8035 _ZN17FractureMechanics10MWLSApprox17getInfluenceNodesEPdb + 1077
>>> [likask:32071] [10] 0   mwls                                0x0000000104f47913 main + 22179
>>> [likask:32071] [11] 0   libdyld.dylib                       0x00007fff904af5ad start + 1
>>> [likask:32071] [12] 0   ???                                 0x000000000000000b 0x0 + 11
>>>
>>>
>>> Can you pleas advise me how I can fix this?
>>>
>>> Kind regards,
>>> Lukasz
>>>
>>> [University of Glasgow: The Times Scottish University of the Year 2018]
>
>