[Gmsh] adaptive meshing: using PView

[walter steffe]

Hello Max, thanks for your replay.

  The mesh field settings follws the example given in gmsh-4.4.1-source/demos/api/viewlist.cpp. The only differencex being that I have used 
scalar tetrahedrons (SS) instead of scalar treiangles (ST). 
>From that example I have understood that the format of data vector has to be:

tet1_V1x, tet1_V2x, tet1_V3x, tet1_V3x, tet1_V4x, tet1_V1y, tet1_V2y, ... tet1_V4z, tet1_val1, tet1_val2, tet1_val3, tet1_val4,
tet2_V1x, ....

Regarding the sharing of my code (which is called EmCAD) it is possible but it would require some work that I was planning to do in the next future.
A few years ago I published an old version of EmCAD (see  https://github.com/wsteffe/EmCAD) but after that I made a lot of changes and I should update it to the
most recent stable version.

The generated pos file is quite big (about 10 MB) but the first two lines are the following:

View "mesh size" {
SS(46.83944890393698,-37.88756912800455,27.48312772362366,47.34769558079194,-37.83372365566942,27.69034188350719,46.76329291980484,-
37.58650617307168,27.60072813488528,46.97847752964434,-37.72705133636931,27.85320246718948){0.6408381837974682};
SS(45.9160848845363,-38.60456203454583,31.84977417926594,46.83102692685463,-38.97962021517143,31.85091832450902,46.00960755639631,-
39.51314340935686,31.79892854402282,46.27960884215373,-38.86612166745061,30.90869736302538){1.042203762895043};


Previously I have tried using a PViewDataGModel based field. 
The printed pos file was the same but the computation of mesh size still wrong.
In that case it failed inside the function PView *getView() const. This function exited after entering inside of the following lines:
  if(v->getData()->hasModel(GModel::current())) {
      Msg::Error(
        "Cannot use view based on current mesh for background mesh: you might"
        " want to use a list-based view (.pos file) instead");
      return 0;
  }

But now I am thinking that probably this was because the current model was not switched to the new one (using gmsh::model::setCurrent())
and was still pointing to the model used for the setting of the PViewDataGModel.

Anyway it is strange that the computation of a field based on the PViewDataList gives wrong results.


Walter



On Sun, 2019-11-17 at 13:06 -0500, Max Orok wrote:
> Hi Walter, 
> 
> My first impression is that it looks like you might indeed be setting the mesh field value using the mesh coordinates. 
> 
> It's a little difficult to debug without an actual runnable program or the pos file that shows everything is OK.
> Would you mind sharing your program and pos file? 
> 
> Sincerely, 
> Max 
> 
> On Sun, Nov 17, 2019 at 4:49 AM walter steffe <walter.steffe at alice.it> wrote:
> > Hello everyone, I am still experimenting with adaptive meshing and I wanted to set a background field based on a PView.
> > I have build the view and related background field with the following code:
> > 
> > 
> >  //sf_ele, and getKeysValues are the same as in gmsh-4.4.1-source/demos/api/adapt_mesh.cpp
> >   getKeysValues(sf_ele, keys, values);
> > 
> > 
> >   int viewTag = gmsh::view::add("mesh size");
> > 
> >   const std::vector<double> listdata;
> >   int nelem=keys.size();
> >   for(int i=0; i<nelem ; i++){
> >     int etag=keys[i];
> >     for(int j=0; j<mesh.elements()[etag]->nodesNum(); j++) listdata.push_back(mesh.elements()[etag]->nodes()[j]->x());
> >     for(int j=0; j<mesh.elements()[etag]->nodesNum(); j++) listdata.push_back(mesh.elements()[etag]->nodes()[j]->y());
> >     for(int j=0; j<mesh.elements()[etag]->nodesNum(); j++) listdata.push_back(mesh.elements()[etag]->nodes()[j]->z());
> >     listdata.insert(listdata.end(),values[i].begin(),values[i].end());
> >   }
> >   gmsh::view::addListData(viewTag, "SS", nelem, listdata);
> > 
> > 
> >   // just to check the data: 
> >   gmsh::view::write(viewTag, "data.pos");
> >   // It seems OK
> > 
> >   ...
> > 
> >   GModel *gm0=new GModel();
> >   #importing of OCC geometry ...
> > 
> >   ...
> > 
> > 
> >   FieldManager *fields = gm->getFields();
> >   int fieldTag=1;
> >   Field *size_f=fields->newField(fieldTag, "PostView");
> >   size_f->options["ViewTag"]->numericalValue(viewTag);
> >   fields->setBackgroundFieldId(fieldTag);
> > 
> > 
> >   gm->mesh(1);     
> >   gm->mesh(2);
> >   gm->mesh(3);
> > 
> > 
> >   The problem is that the field computed in BGM_MeshSize is WRONG.
> > 
> >   Following lines are taken from BackgroundMeshTools.cpp:
> > 
> >   double BGM_MeshSize(GEntity *ge, double U, double V, double X, double Y,
> >                     double Z)
> > {
> >   ....
> >   // lc from fields
> >   double l4 = MAX_LC;
> >   if(ge){
> >     FieldManager *fields = ge->model()->getFields();
> >     if(fields->getBackgroundField() > 0) {
> >       Field *f = fields->get(fields->getBackgroundField());  
> >       if(f) l4 = (*f)(X, Y, Z, ge);
> >     }
> >   }
> > 
> >   ..
> > }
> > 
> > 
> >   I have debugged the code going inside of that computation and I have found that, quite often, it happens that 
> >   the value returned by (*f)(X, Y, Z, ge) coincides with the coordinate (x,y,or z) of a vertex used in the view data.
> > 
> >   This could be produced by a wrong ordering of data passed to gmsh::view::addListData.
> >   But the data file printed by gmsh::view::write(viewTag, "data.pos") seems good and the values are not exchanged with 
> >   the coordinates.
> > 
> >   So I do not understand where is the problem. May you please give me a hint ?
> > 
> > 
> >   Thanks in advance,
> >   Walter Steffè
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> > _______________________________________________
> > gmsh mailing list
> > gmsh at onelab.info
> > http://onelab.info/mailman/listinfo/gmsh
> 
>