Skip to content

Commit

Permalink
Merge pull request #374 from SCOREC/kj/gmsh4
Browse files Browse the repository at this point in the history 
Kj/gmsh4
  • Loading branch information
@cwsmith
cwsmith authored Oct 6, 2022
2 parents 97ca1ef + 7e306ce commit 9992fd4
Show file tree
Hide file tree
Showing 6 changed files with 276 additions and 57 deletions.
8 changes: 7 additions & 1 deletion cmake/FindParmetis.cmake
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,13 @@ if(NOT EXISTS "${METIS_LIBRARY}")
message(FATAL_ERROR "metis library not found")
endif()

set(PARMETIS_LIBRARIES ${PARMETIS_LIBRARY} ${METIS_LIBRARY})
find_library(GK_LIBRARY GKlib PATHS "${PARMETIS_PREFIX}/lib")
if(EXISTS "${GK_LIBRARY}")
set(PARMETIS_LIBRARIES ${PARMETIS_LIBRARY} ${METIS_LIBRARY} ${GK_LIBRARY})
else()
set(PARMETIS_LIBRARIES ${PARMETIS_LIBRARY} ${METIS_LIBRARY})
endif()

set(PARMETIS_INCLUDE_DIRS ${PARMETIS_INCLUDE_DIR} ${METIS_INCLUDE_DIR})

include(FindPackageHandleStandardArgs)
Expand Down
4 changes: 4 additions & 0 deletions mds/apfMDS.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -192,8 +192,12 @@ int getMdsIndex(Mesh2* in, MeshEntity* e);
so call apf::reorderMdsMesh after any mesh modification. */
MeshEntity* getMdsEntity(Mesh2* in, int dimension, int index);

int gmshMajorVersion(const char* filename);

Mesh2* loadMdsFromGmsh(gmi_model* g, const char* filename);

Mesh2* loadMdsDmgFromGmsh(const char* fnameDmg, const char* filename);

Mesh2* loadMdsFromUgrid(gmi_model* g, const char* filename);

void printUgridPtnStats(gmi_model* g, const char* ugridfile, const char* ptnfile,
Expand Down
262 changes: 217 additions & 45 deletions mds/mdsGmsh.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,8 @@ struct Reader {
char* line;
char* word;
size_t linecap;
int major_version;
int minor_version;
bool isQuadratic;
std::map<long, Node> nodeMap;
std::map<long, apf::MeshEntity*> entMap[4];
Expand All @@ -60,19 +62,6 @@ struct Reader {
std::vector<int> physicalType[4];
};

void initReader(Reader* r, apf::Mesh2* m, const char* filename)
{
r->mesh = m;
r->file = fopen(filename, "r");
if (!r->file) {
lion_eprint(1,"couldn't open Gmsh file \"%s\"\n",filename);
abort();
}
r->line = static_cast<char*>(malloc(1));
r->line[0] = '\0';
r->linecap = 1;
r->isQuadratic = false;
}

void freeReader(Reader* r)
{
Expand All @@ -97,6 +86,16 @@ long getLong(Reader* r)
return x;
}

double getDouble(Reader* r)
{
double x;
int nchars;
int ret = sscanf(r->word, "%lf%n", &x, &nchars);
PCU_ALWAYS_ASSERT(ret == 1);
r->word += nchars;
return x;
}

bool startsWith(char const* prefix, char const* s)
{
int ls = strlen(s);
Expand All @@ -118,18 +117,114 @@ void checkMarker(Reader* r, char const* marker)
PCU_ALWAYS_ASSERT(startsWith(marker, r->line));
}

void readNode(Reader* r)
void initReader(Reader* r, apf::Mesh2* m, const char* filename)
{
r->mesh = m;
r->file = fopen(filename, "r");
if (!r->file) {
lion_eprint(1,"couldn't open Gmsh file \"%s\"\n",filename);
abort();
}
r->line = static_cast<char*>(malloc(1));
r->line[0] = '\0';
r->linecap = 1;
r->isQuadratic = false;
seekMarker(r, "$MeshFormat");
int fileType, dataSize;
int ret = sscanf(r->line, "%d.%d %d %d\n",
&r->major_version, &r->minor_version, &fileType, &dataSize);
PCU_ALWAYS_ASSERT(ret==4);
}

void readNode(Reader* r, int bm=-1)
{
long id;
Node n;
apf::Vector3& p = n.point;
sscanf(r->line, "%ld %lf %lf %lf", &id, &p[0], &p[1], &p[2]);
r->nodeMap[id] = n;
if(r->major_version == 4) {
sscanf(r->line, "%lf %lf %lf", &p[0], &p[1], &p[2]);
r->nodeMap[bm] = n;
} else if(r->major_version == 2) {
long id;
sscanf(r->line, "%ld %lf %lf %lf", &id, &p[0], &p[1], &p[2]);
r->nodeMap[id] = n;
}
getLine(r);
}

void readNodes(Reader* r)
void readEntities(Reader* r,const char* fnameDmg)
{
seekMarker(r, "$Entities");
long nlde,ilde,iud,tag,isign,nMV,nME,nMF,nMR;
double x,y,z;
FILE* f = fopen(fnameDmg, "w");
sscanf(r->line, "%ld %ld %ld %ld", &nMV, &nME, &nMF, &nMR);
fprintf(f, "%ld %ld %ld %ld \n", nMR, nMF, nME, nMV); // just reverse order
fprintf(f, "%f %f %f \n ", 0.0, 0.0, 0.0); // Probaby model bounding box?
fprintf(f, "%f %f %f \n", 0.0, 0.0, 0.0); //

getLine(r); // because readNode gets the next line we need this outside for Nodes_Block
for (long i = 0; i < nMV; ++i){
sscanf(r->line, "%ld %lf %lf %lf %ld ", &tag, &x, &y, &z, &iud);
fprintf(f, "%ld %lf %lf %lf \n",tag,x,y,z);
getLine(r);
}
for (long i = 0; i < nME; ++i){
tag = getLong(r);
fprintf(f, "%ld", tag);
for (int i=0; i< 6; ++i) x=getDouble(r); // read past min maxes
iud = getLong(r);
for(long j =0; j < iud; ++j) isign=getLong(r); // read past iud user tags
nlde=getLong(r); // 2 in straight edged models but...
for(long j =0; j < nlde; ++j) {
ilde=getLong(r);
fprintf(f, " %ld", std::abs(ilde)); // modVerts started from 1
}
fprintf(f, "\n");
getLine(r);
}
for (long i = 0; i < nMF; ++i){
tag = getLong(r);
fprintf(f, "%ld %d\n", tag, 1);
for (int i=0; i< 6; ++i) x=getDouble(r); // read past min maxes
iud = getLong(r);
for(long j =0; j < iud; ++j) isign=getLong(r); // read past iud user tags
nlde=getLong(r);
fprintf(f, " %ld \n", nlde);
for(long j =0; j < nlde; ++j) {
ilde=getLong(r);
if(ilde > 0 )
isign=1;
else
isign=0;
fprintf(f, " %ld %ld \n", std::abs(ilde),isign);
}
getLine(r);
}
for (long i = 0; i < nMR; ++i){
tag = getLong(r);
fprintf(f, "%ld %d \n", tag, 1);
for (int i=0; i< 6; ++i) x=getDouble(r); // read past min maxes
iud = getLong(r);
for(long j =0; j < iud; ++j) getLong(r); // read past iud user tags
nlde=getLong(r);
fprintf(f, "%ld \n", nlde);
for(long j =0; j < nlde; ++j) {
ilde=getLong(r);
if(ilde > 0 )
isign=1;
else
isign=0;
fprintf(f, "%ld %ld \n", std::abs(ilde),isign);
}
getLine(r);
}
checkMarker(r, "$EndEntities");
fclose(f);
}

void readNodesV2(Reader* r)
{
PCU_ALWAYS_ASSERT(r->major_version == 2);
seekMarker(r, "$Nodes");
long n = getLong(r);
getLine(r);
Expand All @@ -138,6 +233,27 @@ void readNodes(Reader* r)
checkMarker(r, "$EndNodes");
}

void readNodesV4(Reader* r)
{
PCU_ALWAYS_ASSERT(r->major_version == 4);
seekMarker(r, "$Nodes");
long Num_EntityBlocks,Num_Nodes,Nodes_Block,edim,etag,junk1,junk2,junk3;
sscanf(r->line, "%ld %ld %ld %ld", &Num_EntityBlocks, &Num_Nodes, &junk1, &junk2);
getLine(r); // because readNode gets the next line we need this outside for Nodes_Block
for (long i = 0; i < Num_EntityBlocks; ++i){
sscanf(r->line, "%ld %ld %ld %ld", &edim, &etag, &junk3, &Nodes_Block);
long blockMap[Nodes_Block];
for (long j = 0; j < Nodes_Block; ++j){
getLine(r);
sscanf(r->line, "%ld", &blockMap[j]);
}
getLine(r);
for (long j = 0; j < Nodes_Block; ++j)
readNode(r,blockMap[j]); // has a genLine at end
}
checkMarker(r, "$EndNodes");
}

apf::MeshEntity* lookupVert(Reader* r, long nodeId, apf::ModelEntity* g)
{
PCU_ALWAYS_ASSERT(r->nodeMap.count(nodeId));
Expand All @@ -149,34 +265,38 @@ apf::MeshEntity* lookupVert(Reader* r, long nodeId, apf::ModelEntity* g)
return n.entity;
}

void readElement(Reader* r)
void readElement(Reader* r, long gmshType=-1, long gtag=-1)
{
long id = getLong(r);
long gmshType = getLong(r);
if(r->major_version == 2) {
gmshType = getLong(r);
}
if (isQuadratic(gmshType))
r->isQuadratic = true;
int apfType = apfFromGmsh(gmshType);
PCU_ALWAYS_ASSERT(0 <= apfType);
int nverts = apf::Mesh::adjacentCount[apfType][0];
int dim = apf::Mesh::typeDimension[apfType];
long ntags = getLong(r);
/* The Gmsh 4.9 documentation on the legacy 2.* format states:
* "By default, the first tag is the tag of the physical entity to which the
* element belongs; the second is the tag of the elementary model entity to
* which the element belongs; the third is the number of mesh partitions to
* which the element belongs, followed by the partition ids (negative
* partition ids indicate ghost cells). A zero tag is equivalent to no tag.
* Gmsh and most codes using the MSH 2 format require at least the first two
* tags (physical and elementary tags)."
* A physical entity is a user defined grouping of elementary model entities.
* An elementary model entity is a geometric model entity. */
PCU_ALWAYS_ASSERT(ntags >= 2);
const int physType = static_cast<int>(getLong(r));
PCU_ALWAYS_ASSERT(dim>=0 && dim<4);
r->physicalType[dim].push_back(physType);
long gtag = getLong(r);
for (long i = 2; i < ntags; ++i)
getLong(r); /* discard all other element tags */
if(r->major_version == 2) {
long ntags = getLong(r);
/* The Gmsh 4.9 documentation on the legacy 2.* format states:
* "By default, the first tag is the tag of the physical entity to which the
* element belongs; the second is the tag of the elementary model entity to
* which the element belongs; the third is the number of mesh partitions to
* which the element belongs, followed by the partition ids (negative
* partition ids indicate ghost cells). A zero tag is equivalent to no tag.
* Gmsh and most codes using the MSH 2 format require at least the first two
* tags (physical and elementary tags)."
* A physical entity is a user defined grouping of elementary model entities.
* An elementary model entity is a geometric model entity. */
PCU_ALWAYS_ASSERT(ntags >= 2);
const int physType = static_cast<int>(getLong(r));
PCU_ALWAYS_ASSERT(dim>=0 && dim<4);
r->physicalType[dim].push_back(physType);
gtag = getLong(r);
for (long i = 2; i < ntags; ++i)
getLong(r); /* discard all other element tags */
}
apf::ModelEntity* g = r->mesh->findModelEntity(dim, gtag);
apf::Downward verts;
for (int i = 0; i < nverts; ++i) {
Expand All @@ -193,8 +313,9 @@ void readElement(Reader* r)
getLine(r);
}

void readElements(Reader* r)
void readElementsV2(Reader* r)
{
PCU_ALWAYS_ASSERT(r->major_version == 2);
seekMarker(r, "$Elements");
long n = getLong(r);
getLine(r);
Expand All @@ -203,6 +324,23 @@ void readElements(Reader* r)
checkMarker(r, "$EndElements");
}

void readElementsV4(Reader* r)
{
PCU_ALWAYS_ASSERT(r->major_version == 4);
seekMarker(r, "$Elements");
long Num_EntityBlocks,Num_Elements,Elements_Block,Edim,gtag,gmshType,junk1,junk2;
sscanf(r->line, "%ld %ld %ld %ld", &Num_EntityBlocks, &Num_Elements, &junk1, &junk2);
getLine(r);
for (long i = 0; i < Num_EntityBlocks; ++i){
sscanf(r->line, "%ld %ld %ld %ld", &Edim, &gtag, &gmshType, &Elements_Block);
getLine(r);
for (long j = 0; j < Elements_Block; ++j) {
readElement(r,gmshType,gtag);
}
}
checkMarker(r, "$EndElements");
}

void setElmPhysicalType(Reader* r, apf::Mesh2* m) {
apf::MeshEntity* e;
apf::MeshTag* tag = m->createIntTag("gmsh_physical_entity", 1);
Expand Down Expand Up @@ -281,18 +419,44 @@ void readGmsh(apf::Mesh2* m, const char* filename)
{
Reader r;
initReader(&r, m, filename);
readNodes(&r);
readElements(&r);
m->acceptChanges();
setElmPhysicalType(&r,m);
freeReader(&r);
if(r.major_version == 4) {
readNodesV4(&r);
readElementsV4(&r);
m->acceptChanges();
} else if(r.major_version == 2) {
readNodesV2(&r);
readElementsV2(&r);
m->acceptChanges();
setElmPhysicalType(&r,m);
}
if (r.isQuadratic)
readQuadratic(&r, m, filename);
freeReader(&r);
}
} // closes original namespace

namespace apf {

int gmshMajorVersion(const char* filename) {
Reader r;
Mesh2* m=NULL;
initReader(&r, m, filename);
int version = r.major_version;
freeReader(&r);
return version;
}

void gmshFindDmg(const char* fnameDmg, const char* filename)
{
Reader r;

Mesh2* m=NULL;
initReader(&r, m, filename);
PCU_ALWAYS_ASSERT(r.major_version==4);
readEntities(&r, fnameDmg);
freeReader(&r);
}

namespace apf {

Mesh2* loadMdsFromGmsh(gmi_model* g, const char* filename)
{
Expand All @@ -301,4 +465,12 @@ Mesh2* loadMdsFromGmsh(gmi_model* g, const char* filename)
return m;
}

Mesh2* loadMdsDmgFromGmsh(const char*fnameDmg, const char* filename)
{
gmshFindDmg(fnameDmg, filename); // new function that scans $Entities and writes a dmg
Mesh2* m = makeEmptyMdsMesh(gmi_load(fnameDmg), 0, false);
readGmsh(m, filename);
return m;
}

}
2 changes: 1 addition & 1 deletion pumi-meshes
Submodule pumi-meshes updated 3 files
+88 −0 gmsh/v4/AirfoilDemo.geo
+251,473 −0 gmsh/v4/AirfoilDemo.msh
+151 −0 gmsh/v4/AirfoilDemo_gold.dmg
Loading

0 comments on commit 9992fd4

Please sign in to comment.