[Gmsh] adaptive meshing: using PView
walter steffe
walter.steffe at alice.it
Mon Nov 18 19:04:22 CET 2019
I have found the problem with the PView, it was my fault.
The function getKeysValues(sf_ele, keys, values) I had used is that one defined in gmsh-4.4.1-source/demos/api/adapt_mesh.cpp:
void getKeysValues(const std::map<std::size_t, double> &f,
std::vector<std::size_t> &keys,
std::vector<std::vector<double> > &values)
{
keys.clear();
values.clear();
for(std::map<std::size_t, double>::const_iterator it = f.begin();
it != f.end(); it++){
keys.push_back(it->first);
values.push_back(std::vector<double>(1, it->second));
}
}
The vector associated with a tag (which in my case corresponds with a tetreahedron) has length 1 because the scalar field is constant over the element.
But in my code "listdata.insert(listdata.end(),values[i].begin(),values[i].end())", I was supposing that values[i] contains 4 numbers (the values associated
with the 4 tetrahedron vertices). And my assumption was wrong.
Anyway now I have switched back to the field based on the PViewDataGModel (as done in gmsh-4.4.1-source/demos/api/adapt_mesh.cpp:).
This solution is more efficient (less memory usage) and finally it works well after having set the current model.
Walter
On Mon, 2019-11-18 at 07:13 +0100, walter steffe wrote:
> 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
> >
> >
More information about the gmsh
mailing list