Closed KennethEJansen closed 1 year ago
I guess I found my problem. It wants the 1D elements. changed theEdim==2
to Edim>=1
and it verified.
Since Dan is retired and since Ben drafted me into this SSC challenge, I am drafting him in to help. I am trying to map version 4 gmsh information contained between $Entities and $EndEntities in the .msh file. I did this by hand and got a verified mesh but boundary conditions on model face 0 did not show up for duty in PHASTA which makes me think that one or more places within SCOREC/core there is an assumption that model tags are unique (e.g., you cant have tag=1 on a face AND a vertexz AND a region, AND an edge like gmsh allows. Further, if you are going to use extrusions and the like, gmsh is going to make is own numbering so you can't really control that either.
SO my new function requirements are:
1) Can't load a dmg file because we currently don't have a reliable way to create on from a .msh file (other than Dan's limited cpp route which is orphaned). => gmshUT.cc
(UniqueTags) is forked off of gmsh.sh
/projects/tools/SCOREC-core/core (ExtrusionTagsOffDev)$ ls test/gmsh*.cc
test/gmsh.cc test/gmshUT.cc
23c23,24
< apf::Mesh2* m = apf::loadMdsFromGmsh(gmi_load(argv[1]), argv[2]);
---
> // apf::Mesh2* m = apf::loadMdsFromGmsh(gmi_load(argv[1]), argv[2]);
> apf::Mesh2* m = apf::loadMdsDmgFromGmsh(argv[1], argv[2]);
mds/mdsGmsh.cc
)357 void gmshFindDmg(const char* fnameDmg, const char* filename, int* emap)
358 {
359 Reader r;
360 initReader(&r, m, filename);
361 readEntities(&r, emap);
362 freeReader(r);
363 }
364
365 }
366
367 namespace apf {
368
369 Mesh2* loadMdsFromGmsh(gmi_model* g, const char* filename)
370 {
371 Mesh2* m = makeEmptyMdsMesh(g, 0, false);
372 readGmsh(m, filename);
373 return m;
374 }
375
376 Mesh2* loadMdsDmgFromGmsh(const char*fnameDmg, const char* filename)
377 {
378 int emap[100]; \\letting 100 be the max entities until a stronger C++ coder fixes this
379 gmshFindDmg(fnameDmg, filename, &emap); // new function that scans $Entities and writes a d mg
380 Mesh2* m = makeEmptyMdsMesh(gmi_load(fnameDmg), 0, false);
381 readGmsh(m, filename,*emap);
382 return m;
127 void readEntities(Reader* r,const char* fnameDmg, int* emap, int* ne)
128 {
129 seekMarker(r, "$Entities");
130 long nlde,ilde,iud,junk1,junk2,junk3;
131 double x,y,z
132 FILE* f = fopen(fnameDmg, "w");
133 sscanf(r->line, "%ld %ld %ld %ld", &nme[0], &nme[1], &nme[2], &nme[3]);
134 fprintf(f, "%ld %ld %ld %ld \n", nme[3], nme[2], nme[1], nme[0]); // just reverse order
135 fprintf(f, "%lf %lf %lf \n ", 0, 0, 0); // Probaby model bounding box?
136 fprintf(f, "%lf %lf %lf \n", 0, 0, 0); //
137 getLine(r); // because readNode gets the next line we need this outside for Nodes_Block
138 int entCnt=0;
139 for (long i = 0; i < nme[0]; ++i){
140 sscanf(r->line, "%ld %lf %lf %lf", &tag, &x, &y, &z);
141 emap[entCnt]=tag; // map(dmgTag)=gmshTag is backward how we need but gmshTags dup
142 entCnt++; // vertex entities start from 1
143 fprintf(f, "%ld %lf %lf %lf \n",entCnt,x,y,z);
144 getLine(r);
145 }
146 for (long i = 0; i < nme[1]; ++i){
147 sscanf(r->line, "%ld %lf %lf %lf", &tag, &x, &y, &z);
148 emap[entCnt]=tag;
149 entCnt++;
150 fprintf(f, "%ld", entCnt);
151 sscanf(r->line, "%lf %lf %lf %ld", &x, &y, &z, &iud);
152 for(long j =0; j < iud; j++) sscanf(r->line,ijunk);
153 sscanf(r->line, "%ld", &nlde); // 2 in straight edged models but...
154 for(long j =0; j < nlde; j++) {
155 sscanf(r->line,ilde);
156 for(int k=0; k < nme[0]) { // have to search since map is backwards
157 if(emap[k] == ilde)
158 fprintf(f, "%ld", k+1); // modVerts started from 1
159 }
160 fprintf(f, "\n");
161 getLine(r);
162 }
163 for (long i = 0; i < nme[2]; ++i){
164 sscanf(r->line, "%ld %lf %lf %lf", &tag, &x, &y, &z);
165 emap[entCnt]=tag; // new face tag map
166 entCnt++;
167 fprintf(f, "%ld %ld \n", entCnt, 1);
168 sscanf(r->line, "%lf %lf %lf %ld", &x, &y, &z, &iud);
169 for(long j =0; j < iud; j++) sscanf(r->line,ijunk);
170 sscanf(r->line, "%ld", &nlde);
171 fprintf(f, "%ld \n", nlde);
172 for(long j =0; j < nlde; j++) {
173 sscanf(r->line,ilde);
174 if(ilde > 0 )
175 isign=1;
176 else
177 isign=0;
178 for(int k=nme[0]; k < nme[0]+nme[1]) { // have to search since map is backwards
179 if(emap[k] == ilde)
180 fprintf(f, "%ld %ld \n", k+1,isign);
181 }
182 }
183 getLine(r);184 }
185 for (long i = 0; i < nme[3]; ++i){ //not even sure that this all hangs with nme[3] > 1 but..
186 sscanf(r->line, "%ld %lf %lf %lf", &tag, &x, &y, &z);
187 emap[entCnt]=tag; // new region tag map
188 entCnt++;
189 fprintf(f, "%ld %ld \n", entCnt, 1);
190 sscanf(r->line, "%lf %lf %lf %ld", &x, &y, &z, &iud);
191 for(long j =0; j < iud; j++) sscanf(r->line,ijunk);
192 sscanf(r->line, "%ld", &nlde);
193 fprintf(f, "%ld \n", nlde);
194 for(long j =0; j < nlde; j++) {
195 sscanf(r->line,ilde);
196 if(ilde > 0 )
197 isign=1;
198 else
199 isign=0;
200 for(int k=nme[0]+nme[1]; k < nme[0]+nme[1]+nme[2]) { // have to search since map is back wards
201 if(emap[k] == ilde)
202 fprintf(f, "%ld %ld \n", k+1,isign);
203 }
204 }
205 getLine(r);
206 }
207 checkMarker(r, "$EndEntities");
208 }
Here is the hand made dmg (on the left) with "comments" off to the right of the gmsh .msh file information that the above code meant to read, translate and write out the dmg info (stuff on the left). Of course the new code is also renumbering entities to make them unique and saving the map so that the mds file has a consistent classification (e.g,. run the gmsh data through the map to get mesh classified against unique tags). .dmg .msh (between $Entities and $EndEntities)
1 6 12 8 reverse 8 12 6 1
0 0 0 dmg has two lines with 3 zeros
0 0 0
1 0 -1 0 drop last zero for 8 verts
2 10 -1 0
3 10 1 0
4 0 1 0
5 0 -1 0.1
6 10 -1 0.1
10 10 1 0.1
14 0 1 0.1
1 1 2 1 0 -1 0 10 -1 0 0 2 1 -2
2 2 3 2 10 -1 0 10 1 0 0 2 2 -3
3 3 4 keep keep negate
4 4 1 ends original loop I defined
6 5 6 starts opposite/extruded loop and flows
7 6 10
8 10 14
9 14 5 ends extruded loop
11 1 5 extruded edges in order of root face
12 2 6
16 3 10
20 4 14 last edge
1 1 1 0 -1 0 10 1 0 0 4 1 2 3 4
4 ft skip 7 ne etags in loop order (note this one curls in)
1 1
2 1
3 1
4 1 pos->1
13 1 13 0 -1 0 10 -1 0.1 0 4 1 12 -6 -11
4 neg-> 0
1 1
12 1
6 0
11 0
17 1 17 10 -1 0 10 1 0.1 0 4 2 16 -7 -12
4
2 1
16 1
7 0
12 0
21 1 21 0 1 0 10 1 0.1 0 4 3 20 -8 -16
4
3 1
20 1
8 0
16 0
25 1 25 0 -1 0 0 1 0.1 0 4 4 11 -9 -20
4
4 1
11 1
9 0
20 0
26 1 26 0 -1 0.1 10 1 0.1 0 4 6 7 8 9
4
6 1
7 1
8 1
9 1 nf pos-> 1
1 1 1 0 -1 0 10 1 0.1 0 6 -1 26 13 17 21 25
6 neg-> 0
1 0
26 1
13 1
17 1
21 1
25 1
(base) kjansen@portal1: /projects/tools/SCOREC-core/core (gmshv4Load2MDS)$ git commit -am "made a branch for the not yet finished attempt at a reader from .msh directly into mds making a dmg along the way" [gmshv4Load2MDS cbfc9ce6] made a branch for the not yet finished attempt at a reader from .msh directly into mds making a dmg along the way 2 files changed, 150 insertions(+), 7 deletions(-) create mode 100644 test/gmshUT.cc
In case you would rather track and comment on my efforts here (note however this should be diffed withExtrusionTagsOffDev
Which I don't think has been merged onto develop of SCOREC/core yet even though I did make my changes off of the dev branch at the time (as the name implies).
@cwsmith had mentioned keeping the option to read a V2 gmsh.
That would not be hard for the specific workflow that is described here which builds the dmg and the gmsh .geo file from cpp primitive commands but, this is limited and has no developer so it may be confusing to support that and this direct from .msh approach in the same code. Specifically, v2 does not have the model information in the .msh file that is needed to build the dmg file unless they are built from the same cpp source.
I may not have time to fully dive into this, but it may help to point out that there is Gmsh reader code in Omega_h as well:
https://github.com/sandialabs/omega_h/blob/main/src/Omega_h_gmsh.cpp
It seems to have support for format V4. Perhaps copying this to SCOREC/core would be a path forward?
Soooo the code that I put into PR#372 satisfies verify but based on PHASTAs reaction to the code I am fairly certain that mesh vertices that should be classified on model vertices are not getting that classification.
Further, and possibly the same issue, when I assign a BC to a model face, normal chef protocol is to reconcile through inheritance what should happen on the closure of that model face (e.g., for a model edge it looks at the BCs on both model faces, then for a model vertex, it looks at the BCs on all the model edges incident on that vertex and then through inheritance rules applies a BC). In the model in question, when I do a single layer extrusion of an x-y box in z, the mesh edges at the box corners match model edges (no mesh vertices classified on them) and the mesh vertices (destination and extrusion source) also match model vertices. If I put a comp3 (no slip all three components of velocity) on the faces at the top and bottom of the box, inheritance should apply that same boundary condition to all 8 model vertices each of which should have one mesh vertex classified on it. This is not happening with the committed code. Further it is not happening when bypass inheritance and request in the spj file to put a comp3 on those 8 model vertices directly. Below I will begin to blog my understanding of what the code is doing but I will also pose a more expedient question:
Has anyone tested the original code for its ability to classify mesh on model vertices?
and a second question
Does verify actually check to see if the mesh vertices that are located on model vertices do get classified against them?
If the answer to both is NO, then I think we have to assume the original code is unverified in this capability and I could sure use some help writing some testing code to see if they are correct or not. Or, if it is easier to stare at the code and answer these questions thats fine too.
I might be making a tiny bit of headway understanding what is wrong but I am in so far over my head in terms of understanding the code that does this that the following is likely littered with conceptual miss-understandings. I guess thats what happens when a fortran person tries to use C/C++ code. Hoping someone can see what is going wrong:
gmshFindDmg
which is called from loadMdsDmgFromGmsh
.loadMdsDmgFromGmsh
then calls makeEmptyMdsMesh(gmi_load(fnameDmg),0,false);
that loads the dmg file that gmshFindDmg
just wrote. Since I have visually inspected/checked the dmg file and since that call is the same as the v2.0 version, I have no reason to suspect that I should have needed to change any of that (just moved from test/gmsh.cc down to this level so that is done after the dmg file is created).loadMdsDmgFromGmsh
then calls readGmsh(m,filename,emap);
. Note the new, third argument emap
is a array of integers that stores the gmsh tag for every unique dmg tag. For now I have a hard limit of 100 dmg tags (stored from 0..99) and then use the 4 entries in location 100+d to store how many model d dimensional model entities there are in the model for this mesh. This information allows me to invert the map by doing a while loop search over the appropriate range when I know the gmsh tag and its dim.readGmsh
, I first call initReader
(unchanged)readNodes
. Other than changing the format for how v4 changes the way the data is written and thus must be read there is no change. Specifically, though v4 does have model classification information in each block of nodes it writes, I am not currently using that to classify the mesh while reading the nodes because this was not done in the v2 version of this code and I doubted my ability to do that correctly. Instead, I hoped I could extend what v2 did which was to classify the code during the reading of ELEMENTS (below). Of course I did follow v2's practice of storing the coordinates at n.point and the nodeMap[gmshNode#]=n
which is used in ELEMENTS below. This should be unchanged though it was a little tricky since v4 writes ALL the gmshNode#s classified on a given model entity right after the header that holds the info for that model entity's block of nodes and THEN writes all the coordinate triples after that whereas v2 had that gmshNode# an triple on the same line. I am confident I did this right as the mesh passes verify. readGmsh
then readElements(&r,emap)
. New argument emap is needed to be sure that as I loop over the model entities, each with their own header block that contains
sscanf(r->line, "%ld %ld %ld %ld", &Edim, >ag, &gmshType, &Elements_Block);
dim, gmshtag, type of element, and number of elements in this block, I need emap to map gtag to what I wrote in the dmg file and have currently loaded as described in 6. above. Again, I am confident that this is correct to the level that verify checks.mds/mdsGmsh.cc
)
323 void readElement(Reader* r, long gmshType,long gtag, long skippedLDE)
324 {
325 long id = getLong(r); // tag in 2 and 4
326 id=id-skippedLDE; // id only used for quadratic so this is probably not necessary until that is ready
327 if (isQuadratic(gmshType))
328 r->isQuadratic = true;
329 int apfType = apfFromGmsh(gmshType);
330 PCU_ALWAYS_ASSERT(0 <= apfType);
331 int nverts = apf::Mesh::adjacentCount[apfType][0];
332 int dim = apf::Mesh::typeDimension[apfType];
333 apf::ModelEntity* g = r->mesh->findModelEntity(dim, gtag);
334 apf::Downward verts;
335 for (int i = 0; i < nverts; ++i) {
336 long nid = getLong(r);
337 verts[i] = lookupVert(r, nid, g);
338 }
339 if (dim != 0) {
340 if (dim > r->mesh->getDimension())
341 apf::changeMdsDimension(r->mesh, dim);
342 apf::MeshEntity* ent = apf::buildElement(r->mesh, g, apfType, verts);
343 if (r->isQuadratic)
344 r->entMap[dim][id] = ent;
345 }
346 getLine(r);
347 }
For the first model vertex "element" line 333 does what I would expect. I feed it dim=0 an gtag=1 and it returns g with 1 if I am interpreting TotalView correctly. Great..maybe??. The next model vertex feeds it dim=0 and gtag=2 but now it returns g with a 5. Next gtag=3 returns g with a "9". Then gtag 4 returns g with a d. So it is jumping by 4?? Maybe this is not wrong but I notice that lookupVert(r,nid,g) on line 337 is thinking that this is node 19 if I am interpreting TotalView correctly. Here is a Screenshot of that.
This version, which does process the 0 dim elements but does not add them as elements to the mesh (due to line 339 that was in the original V2 code) does run and pass verify. It runs through Chef as well. PHASTA however does not apply BCs to model vertices which tells me the classification is wrong I think.
Anybody see the problem?
PHASTA seems fine with the current commit for a simple box. However, gmsh dumps construction vertices into the .msh file and when these are not mesh nodes classified on model vertices, the whole premise of verify (and probably SCOREC/core adjacency as a whole) is violated. I confirmed that this can be "fixed" by hand deleting all the spurious entities from mode $Entities, as $Nodes, and as $Elements but this is not a viable solution. I am in touch with Jean Francois Remacle who says the solution is to add Physical Groups. When gmsh sees no Physical Groups, it dumps EVERYTHING. When it sees physical groups, it only dumps the mesh for them.
HOWEVER, while I figure out how to do that I would REALLY appreciate someone to put eyes on readElements and readElement to figure out what that code really needs. Specifically, It obviously needs the highest dimensional elements (volumes in my case) but, can it construct everything else from that or does it also need faces on every model face, edges on every model face and nodes on every model node as I have been giving it (and it has been working with). If so I guess the user will need to make physical groups in the GUI or .geo file for all of these things. If not, I guess the second question is whether the inferred classification that code will create MATCHES the dmg file which is of course required for the downstream workflow.
Another alternative I was considering is to read an auxiliary file that lists the spurious model/mesh vertices that need to be removed and then, at each stage (Entities, Nodes, Elements) check to see if they are "banned" from the mesh and skip over them if they are. This effectively does what I did by deleting them from the .msh file. Tedious and ugly code to be sure. Another solution in the same vein is to scan the elements first to identify the spurious model/mesh vertices (they will be nodes not part of any 1D, 2D, or 3D elements while part of "0D" elements that I have only seen in gmsh).
Yet another solution (might take longer if I don't have help) is to NOT derive model node classification from the $Elements group but instead set it while reading the $Nodes section since it is available there. This would allow us to limit the correct vs. banned check only there as we would only process lines from the $Elements section for dim > 0 and those banned nodes don't appear there. Actually they would still have to not be added to the dmg file in the $Entities section.
Lots of possible solutions but maybe it is good for users to get acquainted with their mesh and set Physical Groups in the GUI or .geo file anyway as some day we can probably replace the spj file with something far more intuitive.
Ok so this is not mainline SCOREC/core stuff and maybe has not been touched in 5 years BUT it would be really cool to keep this capability to load meshes created by gmsh into mds. When I tried to use the tools described on the gmshToPumi page with a recent build of gmsh it failed and so I looked at the gmsh manual to see that this tool was obviously written for a gmsh version 2 mesh which they won't write anymore. Version 4 writes the coordinates in model entity group "blocks" and does the same for elements. I thought I mapped all of READS correctly and it does survive reading the mesh but this crashes on verify for the example case provided at the link above (notched2d is a 2D mesh). Suspects/doubts are: 1) It looks like gmsh starts numbering tags from 1. 2) I am far from certain that I have ported everything that
readElement
is doing inapf
. Here are my code edits.Here is the failed output. Obviously an edge that should be classified on a model edge is classified on a model face but I am not seeing where my code change is goofing that up.