From 09c9b4c3261fd797a63a4df3f7b632413291f28f Mon Sep 17 00:00:00 2001 From: zhang-alvin Date: Fri, 26 Jun 2020 13:30:57 -0400 Subject: [PATCH 01/16] MAINT: remove old comments --- .../MeshAdaptPUMI/createAnalyticGeometry.cpp | 50 +------------------ 1 file changed, 1 insertion(+), 49 deletions(-) diff --git a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp index 1d26aeb4c4..886794d5bf 100644 --- a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp +++ b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp @@ -384,38 +384,9 @@ void makeBox(gmi_model* model) //reparam edges onto face - //int edgeMap[12] = {50,48,46,52,11,16,20,6,73,72,71,74}; int edgeLoop[6][4] = {{50,72,11,73},{72,48,71,16},{46,74,20,71},{52,73,6,74},{50,52,46,48},{11,16,20,6}}; int edgeReparamLoop[6][4] = {{0,3,2,1},{1,0,3,2},{0,1,2,3},{0,1,2,3},{0,1,2,3},{0,3,2,1}}; - //int faceLoop[6] = {80,78,76,82,42,24}; - - //int faceMap[6] = {0,1,2,3,4,5}; - //int faceLoop[6][4] = {{0,9,4,8},{9,1,10,5},{2,11,6,10},{3,8,7,11},{0,3,2,1},{4,5,6,7}}; -/* - typedef void (*ParametricFunctionArray) (double const p[2], double x[3], void*); - ParametricFunctionArray faceFunction[] = - { - face0, - face1, - face2, - face3, - face4, - face5, - }; - for(int i=0; i<6;i++) - { - b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(g_face[i])); - for(int j=0; j<4;j++) - { - agm_use faceUse = add_adj(model, b, edgeLoop[i][j]); - gmi_add_analytic_reparam(model, faceUse, faceFunction[i], 0); - } - } - - -*/ - typedef void (*ParametricFunctionArray) (double const from[2], double to[2], void*); ParametricFunctionArray edgeFaceFunction[] = { @@ -436,14 +407,6 @@ void makeBox(gmi_model* model) } - //make region -/* - int regionPeriodic[3] = {0, 0, 0}; - double regionRanges[3][2] = {{0,0},{0,0},{0.0}}; - gmi_ent* g_region; - g_region = gmi_add_analytic(model, 3, 92, regionFunction, regionPeriodic, regionRanges, 0); - b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(g_region)); -*/ gmi_add_analytic_cell(model,3,92); b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(gmi_find(model,3,92))); @@ -453,18 +416,6 @@ void makeBox(gmi_model* model) gmi_add_analytic_reparam(model, regionUse, regionFunction, 0); } -/* - //add edge adjacencies to region? - it doesn't seem like i need to do this - for(int i=0;i<12;i++) - { - agm_use regionUse = add_adj(model,b,1,edgeMap[i]); - gmi_add_analytic_reparam(model, regionUse, regionFunction, 0); - } -*/ - - //get the sphere reparameterizations onto box - //b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(gmi_find(model,3,92))); agm_use regionUse = add_adj(model, b, 123); gmi_add_analytic_reparam(model, regionUse, regionFunction, 0); @@ -494,6 +445,7 @@ void makeSphere(gmi_model* model) gmi_add_analytic(model, 2, sphereFaceID, sphereFace, faPer, faRan, 0); } +//model tags are based off gmsh default outputs... void setParameterization(gmi_model* model,apf::Mesh2* m) { From 4ea2c391db5501546763d810b2e15282bd7f99e3 Mon Sep 17 00:00:00 2001 From: zhang-alvin Date: Fri, 26 Jun 2020 15:57:08 -0400 Subject: [PATCH 02/16] MAINT: separate adaptation call in analytic geometry creation --- proteus/MeshAdaptPUMI/MeshAdaptPUMI.h | 1 + .../MeshAdaptPUMI/createAnalyticGeometry.cpp | 29 +++++++++++-------- 2 files changed, 18 insertions(+), 12 deletions(-) diff --git a/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h b/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h index ec4c85324b..65382a3461 100644 --- a/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h +++ b/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h @@ -75,6 +75,7 @@ class MeshAdaptPUMIDrvr{ //analytic geometry gmi_model* createSphereInBox(double* boxDim, double*sphereCenter,double radius); void updateSphereCoordinates(double*sphereCenter); + void initialAdapt_analytic(); //Quality Check Functions double getMinimumQuality(); diff --git a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp index 886794d5bf..08f31f6e8a 100644 --- a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp +++ b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp @@ -570,6 +570,9 @@ gmi_model* MeshAdaptPUMIDrvr::createSphereInBox(double* boxDim,double*sphereCent makeSphere(model); + //initial adapt + initialAdapt_analytic(); + //add the box makeBox(model); @@ -577,6 +580,20 @@ gmi_model* MeshAdaptPUMIDrvr::createSphereInBox(double* boxDim,double*sphereCent setParameterization(model,m); m->verify(); + return model; +} + + + +void MeshAdaptPUMIDrvr::updateSphereCoordinates(double*sphereCenter) +{ + xyz_offset[0] = sphereCenter[0]; + xyz_offset[1] = sphereCenter[1]; + xyz_offset[2] = sphereCenter[2]; +} + +void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ + //apf::Field* size_initial = apf::createLagrangeField(m,"size_initial",apf::SCALAR,1); size_iso = apf::createLagrangeField(m,"proteus_size",apf::SCALAR,1); apf::MeshIterator* it = m->begin(0); @@ -633,17 +650,5 @@ gmi_model* MeshAdaptPUMIDrvr::createSphereInBox(double* boxDim,double*sphereCent //apf::writeVtkFiles("initialAdapt2",m); freeField(size_iso); - - return model; } - - -void MeshAdaptPUMIDrvr::updateSphereCoordinates(double*sphereCenter) -{ - xyz_offset[0] = sphereCenter[0]; - xyz_offset[1] = sphereCenter[1]; - xyz_offset[2] = sphereCenter[2]; -} - - From 8c2fe7b8f8bc6d59728fda60ba413bd9e4630432 Mon Sep 17 00:00:00 2001 From: zhang-alvin Date: Mon, 29 Jun 2020 19:16:59 -0400 Subject: [PATCH 03/16] STATE: template for 2D analytic geometry set --- proteus/MeshAdaptPUMI/MeshAdaptPUMI.h | 2 + proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp | 6 +- .../MeshAdaptPUMI/createAnalyticGeometry.cpp | 242 +++++++++++++++--- 3 files changed, 208 insertions(+), 42 deletions(-) diff --git a/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h b/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h index 65382a3461..5e589f130f 100644 --- a/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h +++ b/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h @@ -74,6 +74,8 @@ class MeshAdaptPUMIDrvr{ //analytic geometry gmi_model* createSphereInBox(double* boxDim, double*sphereCenter,double radius); + gmi_model* createCircleInBox(double* boxDim,double*circleCenter, double radius); + void updateSphereCoordinates(double*sphereCenter); void initialAdapt_analytic(); diff --git a/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp b/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp index 0836b538c7..b2ce9280b2 100644 --- a/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp +++ b/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp @@ -164,7 +164,11 @@ int MeshAdaptPUMIDrvr::loadMeshForAnalytic(const char* meshFile,double* boxDim,d m->verify(); //create analytic geometry - gmi_model* testModel = createSphereInBox(boxDim,sphereCenter,radius); + gmi_model* testModel; + if(m->getDimension()==2) + testModel = createCircleInBox(boxDim,sphereCenter,radius); + else + testModel = createSphereInBox(boxDim,sphereCenter,radius); m->verify(); diff --git a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp index 08f31f6e8a..a8a626e36b 100644 --- a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp +++ b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp @@ -423,6 +423,106 @@ void makeBox(gmi_model* model) return; } +class Box{ + public: + void makeBox(gmi_model*); +}; + +void Box::makeBox(gmi_model* model) +{ + //making a box + + int vertPer = 0; + double vertRan[1][2] = {{0.0,0.0}}; + int vertexMap[4] = {58,5,10,56}; + gmi_ent* g_vert[4]; + g_vert[0] = gmi_add_analytic(model, 0, 58, vert0, &vertPer, vertRan, 0); + g_vert[1] = gmi_add_analytic(model, 0, 5, vert1, &vertPer, vertRan, 0); + g_vert[2] = gmi_add_analytic(model, 0, 10, vert2, &vertPer, vertRan, 0); + g_vert[3] = gmi_add_analytic(model, 0, 56, vert3, &vertPer, vertRan, 0); + + int edgePer = 0; + double edgeRan[1][2] = {{0.0,1.0}}; + gmi_ent* g_edge[4]; + + g_edge[0] = gmi_add_analytic(model, 1, 50, edge0, &edgePer, edgeRan, 0); + g_edge[1] = gmi_add_analytic(model, 1, 48, edge1, &edgePer, edgeRan, 0); + g_edge[2] = gmi_add_analytic(model, 1, 46, edge2, &edgePer, edgeRan, 0); + g_edge[3] = gmi_add_analytic(model, 1, 52, edge3, &edgePer, edgeRan, 0); + + //reparameterize vertices on edges + agm_bdry b; + b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(g_edge[0])); + agm_use edgeUse0 = add_adj(model, b, vertexMap[0]); + agm_use edgeUse0_1 = add_adj(model,b,vertexMap[1]); + gmi_add_analytic_reparam(model, edgeUse0, reparamVert_zero, 0); + gmi_add_analytic_reparam(model, edgeUse0_1, reparamVert_one, 0); + + b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(g_edge[1])); + edgeUse0 = add_adj(model, b, vertexMap[1]); + edgeUse0_1 = add_adj(model,b,vertexMap[2]); + gmi_add_analytic_reparam(model, edgeUse0, reparamVert_zero, 0); + gmi_add_analytic_reparam(model, edgeUse0_1, reparamVert_one, 0); + + b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(g_edge[2])); + edgeUse0 = add_adj(model, b, vertexMap[2]); + edgeUse0_1 = add_adj(model,b,vertexMap[3]); + gmi_add_analytic_reparam(model, edgeUse0, reparamVert_one, 0); + gmi_add_analytic_reparam(model, edgeUse0_1, reparamVert_zero, 0); + + b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(g_edge[3])); + edgeUse0 = add_adj(model, b, vertexMap[3]); + edgeUse0_1 = add_adj(model,b,vertexMap[0]); + gmi_add_analytic_reparam(model, edgeUse0, reparamVert_one, 0); + gmi_add_analytic_reparam(model, edgeUse0_1, reparamVert_zero, 0); + + //make faces + + int facePeriodic[2] = {0, 0}; + double faceRanges[2][2] = {{0,1.0},{0,1.0}}; + gmi_ent* g_face[6]; + g_face[0] = gmi_add_analytic(model, 2, 80, face0, facePeriodic, faceRanges, 0); + + //reparam edges onto face + + int edgeLoop[1][4] = {{50,48,46,52}}; + int edgeReparamLoop[1][4] = {{0,1,2,3}}; + + typedef void (*ParametricFunctionArray) (double const from[2], double to[2], void*); + ParametricFunctionArray edgeFaceFunction[] = + { + reparamEdge_0, + reparamEdge_1, + reparamEdge_2, + reparamEdge_3, + }; + + for(int i=0; i<1;i++) + { + b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(g_face[i])); + for(int j=0; j<4;j++) + { + agm_use faceUse = add_adj(model, b, edgeLoop[i][j]); + gmi_add_analytic_reparam(model, faceUse, edgeFaceFunction[edgeReparamLoop[i][j]], 0); + } + } + + + gmi_add_analytic_cell(model,2,92); + + b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(gmi_find(model,2,92))); + for(int i=0; i<1;i++) + { + agm_use regionUse = add_adj(model, b, faceLoop[i]); + gmi_add_analytic_reparam(model, regionUse, regionFunction, 0); + } + + agm_use regionUse = add_adj(model, b, 123); + gmi_add_analytic_reparam(model, regionUse, regionFunction, 0); + + + return; +} //create sphere const double pi = apf::pi; @@ -445,6 +545,33 @@ void makeSphere(gmi_model* model) gmi_add_analytic(model, 2, sphereFaceID, sphereFace, faPer, faRan, 0); } +class Sphere{ + public: + int faceID = 123; + double radius; + double offset[3]; + int dim; + Sphere(int x){ + dim = x; + } + //void sphereFaces(double const*, double*, void*); + //void circleFace(double const, double, void*); + void makeSphere(gmi_model*); + +}; + +void Sphere::makeSphere(gmi_model* model) +{ + int faPer[2] = {1, 0}; + double faRan[2][2] = {{0,6.28318530718},{0.0,apf::pi}}; + if(dim==2){ + faRan[1][1] = 0.0; + } + + sphereRadius = radius; + gmi_add_analytic(model, dim-1, sphereFaceID, sphereFace, faPer, faRan, 0); +} + //model tags are based off gmsh default outputs... void setParameterization(gmi_model* model,apf::Mesh2* m) { @@ -551,47 +678,6 @@ void setParameterization(gmi_model* model,apf::Mesh2* m) m->acceptChanges(); } -gmi_model* MeshAdaptPUMIDrvr::createSphereInBox(double* boxDim,double*sphereCenter, double radius) -{ - sphereRadius = radius; - boxLength = boxDim[0]; - boxWidth = boxDim[1]; - boxHeight = boxDim[2]; - xyz_offset[0] = sphereCenter[0]; - xyz_offset[1] = sphereCenter[1]; - xyz_offset[2] = sphereCenter[2]; - - lion_set_verbosity(1); - - //create the analytic model - gmi_model* model = gmi_make_analytic(); - - //add the sphere - - makeSphere(model); - - //initial adapt - initialAdapt_analytic(); - - //add the box - makeBox(model); - - //apf::writeVtkFiles("initialInitial",m); - setParameterization(model,m); - m->verify(); - - return model; -} - - - -void MeshAdaptPUMIDrvr::updateSphereCoordinates(double*sphereCenter) -{ - xyz_offset[0] = sphereCenter[0]; - xyz_offset[1] = sphereCenter[1]; - xyz_offset[2] = sphereCenter[2]; -} - void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ //apf::Field* size_initial = apf::createLagrangeField(m,"size_initial",apf::SCALAR,1); @@ -652,3 +738,77 @@ void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ freeField(size_iso); } + +void MeshAdaptPUMIDrvr::updateSphereCoordinates(double*sphereCenter) +{ + xyz_offset[0] = sphereCenter[0]; + xyz_offset[1] = sphereCenter[1]; + xyz_offset[2] = sphereCenter[2]; +} + +gmi_model* MeshAdaptPUMIDrvr::createSphereInBox(double* boxDim,double*sphereCenter, double radius) +{ + sphereRadius = radius; + boxLength = boxDim[0]; + boxWidth = boxDim[1]; + boxHeight = boxDim[2]; + xyz_offset[0] = sphereCenter[0]; + xyz_offset[1] = sphereCenter[1]; + xyz_offset[2] = sphereCenter[2]; + + //lion_set_verbosity(1); + + //create the analytic model + gmi_model* model = gmi_make_analytic(); + + //add the sphere + + makeSphere(model); + + //add the box + makeBox(model); + + //initial adapt + initialAdapt_analytic(); + + //apf::writeVtkFiles("initialInitial",m); + setParameterization(model,m); + m->verify(); + + return model; +} + +gmi_model* MeshAdaptPUMIDrvr::createCircleInBox(double* boxDim,double*sphereCenter, double radius) +{ + Sphere* circle = new Sphere(2); + circle->radius = radius; + boxLength = boxDim[0]; + boxWidth = boxDim[1]; + boxHeight = boxDim[2]; + xyz_offset[0] = sphereCenter[0]; + xyz_offset[1] = sphereCenter[1]; + xyz_offset[2] = sphereCenter[2]; + + //lion_set_verbosity(1); + + //create the analytic model + gmi_model* model = gmi_make_analytic(); + + //add the sphere + + circle->makeSphere(model); + + //add the box + Box* box; + box->makeBox(model); + + //initial adapt + initialAdapt_analytic(); + + //apf::writeVtkFiles("initialInitial",m); + setParameterization(model,m); //need to modify + m->verify(); + + return model; +} + From 7c773d8aee3f7255ac8f15f2c05239eee8ffb641 Mon Sep 17 00:00:00 2001 From: zhang-alvin Date: Tue, 30 Jun 2020 19:51:20 -0400 Subject: [PATCH 04/16] STATE: functioning geometry without adapt --- .../MeshAdaptPUMI/createAnalyticGeometry.cpp | 66 +++++++++++++------ 1 file changed, 45 insertions(+), 21 deletions(-) diff --git a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp index a8a626e36b..f40bff7e72 100644 --- a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp +++ b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp @@ -251,9 +251,6 @@ void reparamEdge_3(double const from[2], double to[2], void*) } - - - void regionFunction(double const p[2], double x[3], void*) { (void)p; @@ -497,6 +494,8 @@ void Box::makeBox(gmi_model* model) reparamEdge_3, }; + gmi_add_analytic_cell(model,2,92); + b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(gmi_find(model,2,92))); for(int i=0; i<1;i++) { b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(g_face[i])); @@ -507,16 +506,6 @@ void Box::makeBox(gmi_model* model) } } - - gmi_add_analytic_cell(model,2,92); - - b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(gmi_find(model,2,92))); - for(int i=0; i<1;i++) - { - agm_use regionUse = add_adj(model, b, faceLoop[i]); - gmi_add_analytic_reparam(model, regionUse, regionFunction, 0); - } - agm_use regionUse = add_adj(model, b, 123); gmi_add_analytic_reparam(model, regionUse, regionFunction, 0); @@ -598,6 +587,7 @@ void setParameterization(gmi_model* model,apf::Mesh2* m) } m->setModel(model); m->acceptChanges(); + //Need to set the parametric coordinates of each of the boundary vertices std::map edgeParam; @@ -631,13 +621,13 @@ void setParameterization(gmi_model* model,apf::Mesh2* m) apf::Vector3 newParam; m->getPoint(ent,0,pt); m->getParam(ent,oldParam); - if(modelType==1) + if(modelType==1 && modelTag!=sphereFaceID) { int relevantIndex = edgeParam[modelTag]; newParam[0]=pt[relevantIndex]/edgeLengths[relevantIndex]; m->setParam(ent,newParam); } - else if (modelType==2 && modelTag!=sphereFaceID) + else if (modelType==2 && modelTag!=sphereFaceID && m->getDimension()>2) { int* relevantIndex = faceParam[modelTag][0]; //size is 2 newParam[0] = pt[relevantIndex[0]]/edgeLengths[relevantIndex[0]]; @@ -672,6 +662,36 @@ void setParameterization(gmi_model* model,apf::Mesh2* m) m->setParam(ent,newParam); } + else if (modelType==1 && modelTag == sphereFaceID) //2D version + { + double argy = (pt[1]-xyz_offset[1]); + double argx = (pt[0]-xyz_offset[0]); + if(argx == 0 && argy ==0) + newParam[0] = 0.0; // not sure if this will happen or if this is right + else + newParam[0] = atan2(argy,argx); +/* + double arg2 = (pt[2]-xyz_offset[2])/sphereRadius; + if(arg2 < -1.0) + arg2 = -1.0; + else if (arg2 > 1.0) + arg2 = 1.0; + + newParam[1] = acos(arg2); + if(newParam[0]<0) + newParam[0] = newParam[0]+2*apf::pi; + if(newParam[0]>2*apf::pi) + newParam[0] = newParam[0]-2*apf::pi; + //this is probably unnecessary + if(newParam[1]<0.0) + newParam[1] = -1*newParam[1]; + if(newParam[1]>apf::pi) + newParam[1] = -1*(newParam[1]-2.0*apf::pi); +*/ + + m->setParam(ent,newParam); + } + } //end if } //end while m->end(it); @@ -697,7 +717,7 @@ void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ gradeMesh(1.5); - + std::cout<<"adapt0\n"; ma::Input* in = ma::configure(m,size_iso); in->maximumIterations = 10; in->shouldSnap = true; @@ -705,6 +725,7 @@ void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ in->shouldFixShape = true; in->debugFolder="./debug_fine"; ma::adaptVerbose(in,false); + std::cout<<"adapt1\n"; m->verify(); //apf::writeVtkFiles("initialAdapt",m); @@ -768,13 +789,14 @@ gmi_model* MeshAdaptPUMIDrvr::createSphereInBox(double* boxDim,double*sphereCent //add the box makeBox(model); - //initial adapt - initialAdapt_analytic(); - //apf::writeVtkFiles("initialInitial",m); setParameterization(model,m); m->verify(); + //initial adapt + initialAdapt_analytic(); + + m->verify(); return model; } @@ -802,11 +824,13 @@ gmi_model* MeshAdaptPUMIDrvr::createCircleInBox(double* boxDim,double*sphereCent Box* box; box->makeBox(model); + //apf::writeVtkFiles("initialInitial",m); + setParameterization(model,m); //need to modify + m->verify(); + //initial adapt initialAdapt_analytic(); - //apf::writeVtkFiles("initialInitial",m); - setParameterization(model,m); //need to modify m->verify(); return model; From 6a9c5146d437678040a31041e900a1df3288e204 Mon Sep 17 00:00:00 2001 From: zhang-alvin Date: Thu, 2 Jul 2020 18:37:20 -0400 Subject: [PATCH 05/16] STATE: functioning box-only config --- .../MeshAdaptPUMI/createAnalyticGeometry.cpp | 169 +++++++++++++----- 1 file changed, 120 insertions(+), 49 deletions(-) diff --git a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp index f40bff7e72..e8c47ebbd2 100644 --- a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp +++ b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp @@ -236,17 +236,17 @@ void reparamEdge_0(double const from[2], double to[2], void*) void reparamEdge_1(double const from[2], double to[2], void*) { to[0] = 0.0; - to[1] = from[0]; + to[1] = 1.0-from[0]; } void reparamEdge_2(double const from[2], double to[2], void*) { - to[0] = from[0]; + to[0] = 1.0 - from[0]; to[1] = 1.0; } void reparamEdge_3(double const from[2], double to[2], void*) { - to[0] = 1.0; + to[0] = 1.0; to[1] = from[0]; } @@ -464,52 +464,42 @@ void Box::makeBox(gmi_model* model) b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(g_edge[2])); edgeUse0 = add_adj(model, b, vertexMap[2]); edgeUse0_1 = add_adj(model,b,vertexMap[3]); - gmi_add_analytic_reparam(model, edgeUse0, reparamVert_one, 0); - gmi_add_analytic_reparam(model, edgeUse0_1, reparamVert_zero, 0); + gmi_add_analytic_reparam(model, edgeUse0, reparamVert_zero, 0); + gmi_add_analytic_reparam(model, edgeUse0_1, reparamVert_one, 0); b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(g_edge[3])); edgeUse0 = add_adj(model, b, vertexMap[3]); edgeUse0_1 = add_adj(model,b,vertexMap[0]); - gmi_add_analytic_reparam(model, edgeUse0, reparamVert_one, 0); - gmi_add_analytic_reparam(model, edgeUse0_1, reparamVert_zero, 0); - - //make faces - - int facePeriodic[2] = {0, 0}; - double faceRanges[2][2] = {{0,1.0},{0,1.0}}; - gmi_ent* g_face[6]; - g_face[0] = gmi_add_analytic(model, 2, 80, face0, facePeriodic, faceRanges, 0); + gmi_add_analytic_reparam(model, edgeUse0, reparamVert_zero, 0); + gmi_add_analytic_reparam(model, edgeUse0_1, reparamVert_one, 0); //reparam edges onto face + int edgeLoop[4] = {50,48,46,52}; + + //gmi_add_analytic_cell(model,2,92); + int faPer[2] = {0, 0}; + double faRan[2][2] = {{0,1},{0,1}}; - int edgeLoop[1][4] = {{50,48,46,52}}; - int edgeReparamLoop[1][4] = {{0,1,2,3}}; - typedef void (*ParametricFunctionArray) (double const from[2], double to[2], void*); ParametricFunctionArray edgeFaceFunction[] = { reparamEdge_0, - reparamEdge_1, - reparamEdge_2, reparamEdge_3, + reparamEdge_2, + reparamEdge_1, }; - gmi_add_analytic_cell(model,2,92); + gmi_ent* f = gmi_add_analytic(model, 2, 92, face4, faPer, faRan, 0); b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(gmi_find(model,2,92))); - for(int i=0; i<1;i++) + for(int i=0; i<4;i++) { - b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(g_face[i])); - for(int j=0; j<4;j++) - { - agm_use faceUse = add_adj(model, b, edgeLoop[i][j]); - gmi_add_analytic_reparam(model, faceUse, edgeFaceFunction[edgeReparamLoop[i][j]], 0); - } + agm_use regionUse = add_adj(model, b, edgeLoop[i]); + gmi_add_analytic_reparam(model, regionUse, edgeFaceFunction[i], 0); } agm_use regionUse = add_adj(model, b, 123); gmi_add_analytic_reparam(model, regionUse, regionFunction, 0); - return; } @@ -670,27 +660,94 @@ void setParameterization(gmi_model* model,apf::Mesh2* m) newParam[0] = 0.0; // not sure if this will happen or if this is right else newParam[0] = atan2(argy,argx); -/* - double arg2 = (pt[2]-xyz_offset[2])/sphereRadius; - if(arg2 < -1.0) - arg2 = -1.0; - else if (arg2 > 1.0) - arg2 = 1.0; - newParam[1] = acos(arg2); - if(newParam[0]<0) - newParam[0] = newParam[0]+2*apf::pi; - if(newParam[0]>2*apf::pi) - newParam[0] = newParam[0]-2*apf::pi; - //this is probably unnecessary - if(newParam[1]<0.0) - newParam[1] = -1*newParam[1]; - if(newParam[1]>apf::pi) - newParam[1] = -1*(newParam[1]-2.0*apf::pi); -*/ + m->setParam(ent,newParam); + } + + } //end if + } //end while + m->end(it); + m->acceptChanges(); +} + +void setParameterization2D(gmi_model* model,apf::Mesh2* m) +{ + //Get the classification of each entity in the SimMesh + apf::MeshEntity* ent; + for(int i =0;i<3;i++) + { + apf::MeshIterator* it = m->begin(i); + while( (ent = m->iterate(it))) + { + apf::ModelEntity* g_ent = m->toModel(ent); + int modelTag = m->getModelTag(g_ent); + int modelType = m->getModelType(g_ent); + if(modelType==2 && modelTag!=92) + std::cout<<"model tag,type "< 139 && modelTag< 148) + { + m->setModelEntity(ent,(apf::ModelEntity*)gmi_find(model,modelType,sphereFaceID)); + } + else + m->setModelEntity(ent,(apf::ModelEntity*)gmi_find(model,modelType,modelTag)); + + } + m->end(it); + } + m->setModel(model); + m->acceptChanges(); + + //Need to set the parametric coordinates of each of the boundary vertices + std::map edgeParam; + const int numEdges = 4; + int edgeScales[numEdges] = {0,1,0,1}; + double edgeLengths[2] = {boxLength,boxWidth}; + for(int i=0;ibegin(0); + while( (ent = m->iterate(it)) ) + { + apf::ModelEntity* g_ent = m->toModel(ent); + + apf::MeshEntity* ev[2]; + m->getDownward(ent,0,ev); + int modelTag = m->getModelTag(g_ent); + int modelType = m->getModelType(g_ent); + if(modelType<3 && modelType!=0) + { + apf::Vector3 pt; + apf::Vector3 oldParam; + apf::Vector3 newParam(0.0,0.0,0.0); + m->getPoint(ent,0,pt); + m->getParam(ent,oldParam); + if(modelType==1 && modelTag!=sphereFaceID) + { + int relevantIndex = edgeParam[modelTag]; + newParam[0]=pt[relevantIndex]/edgeLengths[relevantIndex]; + std::cout<<"model tag "<setParam(ent,newParam); + } + else if (modelType==1 && modelTag == sphereFaceID) //2D version + { + std::cout<<"xyz offset "<setParam(ent,newParam); } + else if (modelType==2) + { + newParam[0] = pt[0]/boxLength; + newParam[1] = pt[1]/boxWidth; + m->setParam(ent,newParam); + } } //end if } //end while @@ -704,20 +761,32 @@ void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ size_iso = apf::createLagrangeField(m,"proteus_size",apf::SCALAR,1); apf::MeshIterator* it = m->begin(0); apf::MeshEntity* ent; + hmin = 0.04; + hmax = 0.1; + std::cout<<"hmin, hmax "<iterate(it)) ) { +/* apf::Vector3 pt; m->getPoint(ent,0,pt); if(sqrt( (pt[0]-xyz_offset[0])*(pt[0]-xyz_offset[0])+ (pt[1]-xyz_offset[1])*(pt[1]-xyz_offset[1]) + (pt[2]-xyz_offset[2])*(pt[2]-xyz_offset[2])) < sphereRadius*1.5) + { apf::setScalar(size_iso,ent,0,hmin); + } else + { apf::setScalar(size_iso,ent,0,hmax); + } +*/ + apf::setScalar(size_iso,ent,0,hmin); } m->end(it); - gradeMesh(1.5); + //gradeMesh(1.5); std::cout<<"adapt0\n"; + apf::writeVtkFiles("initialProteus",m); ma::Input* in = ma::configure(m,size_iso); in->maximumIterations = 10; in->shouldSnap = true; @@ -731,6 +800,7 @@ void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ //apf::writeVtkFiles("initialAdapt",m); freeField(size_iso); +/* size_iso = apf::createLagrangeField(m,"proteus_size",apf::SCALAR,1); it = m->begin(0); while( (ent = m->iterate(it)) ) @@ -757,6 +827,7 @@ void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ //apf::writeVtkFiles("initialAdapt2",m); freeField(size_iso); +*/ } @@ -811,21 +882,21 @@ gmi_model* MeshAdaptPUMIDrvr::createCircleInBox(double* boxDim,double*sphereCent xyz_offset[1] = sphereCenter[1]; xyz_offset[2] = sphereCenter[2]; - //lion_set_verbosity(1); + lion_set_verbosity(1); //create the analytic model gmi_model* model = gmi_make_analytic(); //add the sphere - circle->makeSphere(model); + //circle->makeSphere(model); //add the box Box* box; box->makeBox(model); //apf::writeVtkFiles("initialInitial",m); - setParameterization(model,m); //need to modify + setParameterization2D(model,m); //need to modify m->verify(); //initial adapt From d6b1366e27b8432bc5f5343af82d7704be7decc0 Mon Sep 17 00:00:00 2001 From: zhang-alvin Date: Mon, 6 Jul 2020 13:42:56 -0400 Subject: [PATCH 06/16] STATE: functioning circle --- .../MeshAdaptPUMI/createAnalyticGeometry.cpp | 65 ++++++++++++++++--- 1 file changed, 55 insertions(+), 10 deletions(-) diff --git a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp index e8c47ebbd2..1033e11f1f 100644 --- a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp +++ b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp @@ -257,6 +257,44 @@ void regionFunction(double const p[2], double x[3], void*) (void)x; } +//from edge parameterization to face parameterization +/* +void reparam_Circle(double const from[2], double to[2], void*){ + double x = from[0]*boxLength; + double y = from[1]*boxWidth; + std::cout<<"from[0],from[1],x,y "<begin(0); apf::MeshEntity* ent; - hmin = 0.04; - hmax = 0.1; - std::cout<<"hmin, hmax "<iterate(it)) ) { -/* apf::Vector3 pt; m->getPoint(ent,0,pt); if(sqrt( (pt[0]-xyz_offset[0])*(pt[0]-xyz_offset[0])+ (pt[1]-xyz_offset[1])*(pt[1]-xyz_offset[1]) + (pt[2]-xyz_offset[2])*(pt[2]-xyz_offset[2])) < sphereRadius*1.5) @@ -778,8 +824,6 @@ void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ { apf::setScalar(size_iso,ent,0,hmax); } -*/ - apf::setScalar(size_iso,ent,0,hmin); } m->end(it); @@ -796,6 +840,7 @@ void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ ma::adaptVerbose(in,false); std::cout<<"adapt1\n"; m->verify(); + apf::writeVtkFiles("finalProteus",m); //apf::writeVtkFiles("initialAdapt",m); freeField(size_iso); @@ -889,7 +934,7 @@ gmi_model* MeshAdaptPUMIDrvr::createCircleInBox(double* boxDim,double*sphereCent //add the sphere - //circle->makeSphere(model); + circle->makeSphere(model); //add the box Box* box; From 45bacacc9b0d3f1c8ed89929ae32103e9d0ee1df Mon Sep 17 00:00:00 2001 From: zhang-alvin Date: Mon, 6 Jul 2020 15:47:56 -0400 Subject: [PATCH 07/16] ENH: add isAnalytic flag to determine whether snapping is used --- proteus/MeshAdaptPUMI/MeshAdaptPUMI.h | 1 + proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp | 4 ++++ proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp | 15 +++++++-------- 3 files changed, 12 insertions(+), 8 deletions(-) diff --git a/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h b/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h index 5e589f130f..86b6e7b874 100644 --- a/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h +++ b/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h @@ -110,6 +110,7 @@ class MeshAdaptPUMIDrvr{ bool hasAnalyticSphere; bool useProteus; bool useProteusAniso; + bool isAnalytic; diff --git a/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp b/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp index b2ce9280b2..32b9e0a255 100644 --- a/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp +++ b/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp @@ -84,6 +84,7 @@ MeshAdaptPUMIDrvr::MeshAdaptPUMIDrvr() initialReconstructed = 0; maxAspect = 2.0;//maxAspectRatio; hasIBM=hasInterface=hasVMS=hasERM=hasAniso=hasAnalyticSphere=useProteus=useProteusAniso=0; + isAnalytic=0; } MeshAdaptPUMIDrvr::~MeshAdaptPUMIDrvr() @@ -797,6 +798,7 @@ int MeshAdaptPUMIDrvr::adaptPUMIMesh(const char* inputString) in->shouldRunPostZoltan = true; in->maximumImbalance = 1.05; in->maximumIterations = numIter; +/* if(size_field_config == "meshQuality") { in->shouldSnap = true; @@ -804,6 +806,8 @@ int MeshAdaptPUMIDrvr::adaptPUMIMesh(const char* inputString) } else in->shouldSnap = false; +*/ + in->shouldSnap=isAnalytic; //in->goodQuality = 0.16;//0.027; //double mass_before = getTotalMass(); diff --git a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp index 1033e11f1f..0ed5fd96e5 100644 --- a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp +++ b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp @@ -289,8 +289,8 @@ void reparam_Circle(double const from[2], double to[2], void*){ double y = sphereRadius*sin(from[0])+xyz_offset[1]; to[0] = x/boxLength; to[1] = y/boxWidth; - std::cout<<"box length width "<< boxLength<<" "<toModel(ent); int modelTag = m->getModelTag(g_ent); int modelType = m->getModelType(g_ent); - if(modelType==2 && modelTag!=92) - std::cout<<"model tag,type "< 139 && modelTag< 148) { m->setModelEntity(ent,(apf::ModelEntity*)gmi_find(model,modelType,sphereFaceID)); @@ -778,12 +776,12 @@ void setParameterization2D(gmi_model* model,apf::Mesh2* m) { int relevantIndex = edgeParam[modelTag]; newParam[0]=pt[relevantIndex]/edgeLengths[relevantIndex]; - std::cout<<"model tag "<setParam(ent,newParam); } else if (modelType==1 && modelTag == sphereFaceID) //2D version { - std::cout<<"xyz offset "<begin(0); @@ -945,7 +943,8 @@ gmi_model* MeshAdaptPUMIDrvr::createCircleInBox(double* boxDim,double*sphereCent m->verify(); //initial adapt - initialAdapt_analytic(); + //initialAdapt_analytic(); + isAnalytic=1; m->verify(); From aec29b10d79448c4d8c4b281564db83c81a3a8bd Mon Sep 17 00:00:00 2001 From: zhang-alvin Date: Tue, 14 Jul 2020 12:21:22 -0400 Subject: [PATCH 08/16] ENH: allow for users to set barycenters, which is needed for adaptivity hack --- proteus/TwoPhaseFlow/utils/Parameters.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/proteus/TwoPhaseFlow/utils/Parameters.py b/proteus/TwoPhaseFlow/utils/Parameters.py index 64aa9abc07..32bbf21c92 100644 --- a/proteus/TwoPhaseFlow/utils/Parameters.py +++ b/proteus/TwoPhaseFlow/utils/Parameters.py @@ -1730,7 +1730,8 @@ def _initializePhysics(self): # COEFFICIENTS coeffs = self.p.coefficients coeffs.flowModelIndex = V_model - coeffs.barycenters = domain.barycenters + if(coeffs.barycenters is None): + coeffs.barycenters = domain.barycenters coeffs.nd = nd coeffs.initialize() # INITIAL CONDITIONS From da307e747fc17455fad2241744604f55e9d9387a Mon Sep 17 00:00:00 2001 From: zhang-alvin Date: Tue, 14 Jul 2020 13:52:39 -0400 Subject: [PATCH 09/16] STATE: before refactoring into factory pattern --- .../MeshAdaptPUMI/createAnalyticGeometry.cpp | 30 +++++++++---------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp index 0ed5fd96e5..664ee378d3 100644 --- a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp +++ b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp @@ -293,9 +293,6 @@ void reparam_Circle(double const from[2], double to[2], void*){ //std::cout<<"from "<toModel(ent); int modelTag = m->getModelTag(g_ent); int modelType = m->getModelType(g_ent); + //std::cout<<"modelTag "< 139 && modelTag< 148) { m->setModelEntity(ent,(apf::ModelEntity*)gmi_find(model,modelType,sphereFaceID)); @@ -814,6 +812,7 @@ void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ { apf::Vector3 pt; m->getPoint(ent,0,pt); +/* if(sqrt( (pt[0]-xyz_offset[0])*(pt[0]-xyz_offset[0])+ (pt[1]-xyz_offset[1])*(pt[1]-xyz_offset[1]) + (pt[2]-xyz_offset[2])*(pt[2]-xyz_offset[2])) < sphereRadius*1.5) { apf::setScalar(size_iso,ent,0,hmin); @@ -822,12 +821,13 @@ void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ { apf::setScalar(size_iso,ent,0,hmax); } +*/ + apf::setScalar(size_iso,ent,0,0.1); } m->end(it); //gradeMesh(1.5); - std::cout<<"adapt0\n"; apf::writeVtkFiles("initialProteus",m); ma::Input* in = ma::configure(m,size_iso); in->maximumIterations = 10; @@ -836,28 +836,29 @@ void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ in->shouldFixShape = true; in->debugFolder="./debug_fine"; ma::adaptVerbose(in,false); - std::cout<<"adapt1\n"; m->verify(); - apf::writeVtkFiles("finalProteus",m); + apf::writeVtkFiles("middleProteus",m); //apf::writeVtkFiles("initialAdapt",m); freeField(size_iso); -/* size_iso = apf::createLagrangeField(m,"proteus_size",apf::SCALAR,1); it = m->begin(0); while( (ent = m->iterate(it)) ) { +/* apf::Vector3 pt; m->getPoint(ent,0,pt); if(sqrt( (pt[0]-xyz_offset[0])*(pt[0]-xyz_offset[0])+ (pt[1]-xyz_offset[1])*(pt[1]-xyz_offset[1]) + (pt[2]-xyz_offset[2])*(pt[2]-xyz_offset[2])) < sphereRadius*1.5) apf::setScalar(size_iso,ent,0,hmin); else apf::setScalar(size_iso,ent,0,hmax); +*/ + apf::setScalar(size_iso,ent,0,0.01); } m->end(it); - gradeMesh(1.5); + //gradeMesh(1.5); in = ma::configure(m,size_iso); in->maximumIterations = 10; @@ -868,9 +869,9 @@ void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ ma::adaptVerbose(in,false); m->verify(); + apf::writeVtkFiles("finalProteus",m); //apf::writeVtkFiles("initialAdapt2",m); freeField(size_iso); -*/ } @@ -916,8 +917,8 @@ gmi_model* MeshAdaptPUMIDrvr::createSphereInBox(double* boxDim,double*sphereCent gmi_model* MeshAdaptPUMIDrvr::createCircleInBox(double* boxDim,double*sphereCenter, double radius) { - Sphere* circle = new Sphere(2); - circle->radius = radius; + Sphere circle = Sphere(2); + circle.radius = radius; boxLength = boxDim[0]; boxWidth = boxDim[1]; boxHeight = boxDim[2]; @@ -932,13 +933,12 @@ gmi_model* MeshAdaptPUMIDrvr::createCircleInBox(double* boxDim,double*sphereCent //add the sphere - circle->makeSphere(model); + circle.makeSphere(model); //add the box - Box* box; - box->makeBox(model); + Box box; + box.makeBox(model); - //apf::writeVtkFiles("initialInitial",m); setParameterization2D(model,m); //need to modify m->verify(); From e58bf56bbddc81750989c49f9274a822c8382850 Mon Sep 17 00:00:00 2001 From: zhang-alvin Date: Tue, 14 Jul 2020 16:15:21 -0400 Subject: [PATCH 10/16] MAINT: simplified 2D box creation --- .../MeshAdaptPUMI/createAnalyticGeometry.cpp | 153 ++++++++---------- 1 file changed, 69 insertions(+), 84 deletions(-) diff --git a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp index 664ee378d3..e965a089f6 100644 --- a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp +++ b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp @@ -257,31 +257,7 @@ void regionFunction(double const p[2], double x[3], void*) (void)x; } -//from edge parameterization to face parameterization -/* -void reparam_Circle(double const from[2], double to[2], void*){ - double x = from[0]*boxLength; - double y = from[1]*boxWidth; - std::cout<<"from[0],from[1],x,y "< vertexMap; + double vertRan[1][2]={{0.0,0.0}}; + int vertPer=0; + + std::vector edgeMap; + int edgePer = 0; + double edgeRan[1][2] = {{0.0,1.0}}; + + int faPer[2] = {0, 0}; + double faRan[2][2] = {{0,1},{0,1}}; + int regionID = 92; // fixed ID +}; + +/* +class AnalyticGeom(int dim, info){ + + //Box; + //Sphere; + public: + model createGeometry() + create2Dbox; + create3Dbox; +*/ +class Box : public Enclosure{ public: void makeBox(gmi_model*); }; @@ -463,59 +465,45 @@ class Box{ void Box::makeBox(gmi_model* model) { //making a box + vertexMap = {58,5,10,56}; + gmi_ent* g_vert[vertexMap.size()]; + + EntityMapArray vertexPoints[] = { + vert0, + vert1, + vert2, + vert3 + }; - int vertPer = 0; - double vertRan[1][2] = {{0.0,0.0}}; - int vertexMap[4] = {58,5,10,56}; - gmi_ent* g_vert[4]; - g_vert[0] = gmi_add_analytic(model, 0, 58, vert0, &vertPer, vertRan, 0); - g_vert[1] = gmi_add_analytic(model, 0, 5, vert1, &vertPer, vertRan, 0); - g_vert[2] = gmi_add_analytic(model, 0, 10, vert2, &vertPer, vertRan, 0); - g_vert[3] = gmi_add_analytic(model, 0, 56, vert3, &vertPer, vertRan, 0); + for(auto i=0; i> indexPairs = {{0,1}, {1,2}, {2,3}, {3,0}}; - //reparam edges onto face - int edgeLoop[4] = {50,48,46,52}; + for(int i=0;iverify(); From 4d522f5823109cd585481bd9e4ce749cc45a9405 Mon Sep 17 00:00:00 2001 From: zhang-alvin Date: Wed, 15 Jul 2020 15:17:29 -0400 Subject: [PATCH 11/16] MAINT: cleaned up analytic geometry creation and generalized structures for cleaner code. unclear if 3D functionality works --- proteus/MeshAdaptPUMI/MeshAdaptPUMI.h | 3 +- proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp | 16 +- .../MeshAdaptPUMI/createAnalyticGeometry.cpp | 886 ++++++++---------- .../MeshAdaptPUMI/createAnalyticGeometry.h | 13 - 4 files changed, 393 insertions(+), 525 deletions(-) diff --git a/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h b/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h index 86b6e7b874..2180558a5f 100644 --- a/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h +++ b/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h @@ -73,8 +73,7 @@ class MeshAdaptPUMIDrvr{ int gradeMesh(double gradationFactor); //analytic geometry - gmi_model* createSphereInBox(double* boxDim, double*sphereCenter,double radius); - gmi_model* createCircleInBox(double* boxDim,double*circleCenter, double radius); + void createAnalyticGeometry(int dim, double* boxDim, double*sphereCenter,double radius); void updateSphereCoordinates(double*sphereCenter); void initialAdapt_analytic(); diff --git a/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp b/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp index 32b9e0a255..3ceaec3695 100644 --- a/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp +++ b/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp @@ -161,23 +161,13 @@ int MeshAdaptPUMIDrvr::loadMeshForAnalytic(const char* meshFile,double* boxDim,d //assume analytic comm_size = PCU_Comm_Peers(); comm_rank = PCU_Comm_Self(); - m = apf::loadMdsMesh(".null", meshFile); + //m = apf::loadMdsMesh(".null", meshFile); + m = apf::loadMdsFromGmsh(gmi_load(".null"), meshFile); m->verify(); //create analytic geometry - gmi_model* testModel; - if(m->getDimension()==2) - testModel = createCircleInBox(boxDim,sphereCenter,radius); - else - testModel = createSphereInBox(boxDim,sphereCenter,radius); + createAnalyticGeometry(m->getDimension(),boxDim,sphereCenter,radius); m->verify(); - - -/* - apf::writeVtkFiles("afterAnalytic",m); - std::cout<<"test Model "<getModel()< faceMap; int faPer[2] = {0, 0}; double faRan[2][2] = {{0,1},{0,1}}; - int regionID = 92; // fixed ID -}; - -/* -class AnalyticGeom(int dim, info){ - //Box; - //Sphere; - public: - model createGeometry() - create2Dbox; - create3Dbox; -*/ -class Box : public Enclosure{ - public: - void makeBox(gmi_model*); + int regionID = 92; // fixed ID + void makeBox2D(gmi_model* model); + void makeBox3D(gmi_model* model); }; -void Box::makeBox(gmi_model* model) +void Enclosure::makeBox2D(gmi_model* model) { //making a box vertexMap = {58,5,10,56}; @@ -523,64 +378,144 @@ void Box::makeBox(gmi_model* model) return; } -//create sphere -const double pi = apf::pi; -int sphereFaceID = 123; -double sphereRadius; -double xyz_offset[3]; - -void sphereFace(double const p[2], double x[3], void*) +void Enclosure::makeBox3D(gmi_model* model) { - x[0] = xyz_offset[0]+sphereRadius*cos(p[0]) * sin(p[1]); - x[1] = xyz_offset[1]+sphereRadius*sin(p[0]) * sin(p[1]); - x[2] = xyz_offset[2]+sphereRadius*cos(p[1]); -} + //making a box -void sphereFace2D(double const p[2], double x[3], void*) -{ - x[0] = xyz_offset[0]+sphereRadius*cos(p[0]); - x[1] = xyz_offset[1]+sphereRadius*sin(p[0]); - x[2] = 0.0; -} + vertexMap = {58,56,54,60,5,10,15,2}; + gmi_ent* g_vert[vertexMap.size()]; + EntityMapArray vertexPoints[] = { + vert0, + vert1, + vert2, + vert3, + vert4, + vert5, + vert6, + vert7 + }; -void makeSphere(gmi_model* model) -{ - int faPer[2] = {1, 0}; - double faRan[2][2] = {{0,6.28318530718},{0.0,apf::pi}}; + for(auto i=0; i> indexPairs = {{0,1}, {1,2}, {2,3}, {3,0}, + {4,5},{5,6},{6,7},{7,4},{0,4},{1,5},{2,6},{3,7} }; + + for(int i=0;igetModelType(g_ent); if(modelTag > 139 && modelTag< 148) { - m->setModelEntity(ent,(apf::ModelEntity*)gmi_find(model,modelType,sphereFaceID)); + m->setModelEntity(ent,(apf::ModelEntity*)gmi_find(model,modelType,sphere.faceID)); } else m->setModelEntity(ent,(apf::ModelEntity*)gmi_find(model,modelType,modelTag)); @@ -613,14 +548,14 @@ void setParameterization(gmi_model* model,apf::Mesh2* m) double edgeLengths[3] = {boxLength,boxWidth,boxHeight}; for(int i=0;i<12;i++) { - edgeParam[edgeMap[i]] = edgeScales[i]; + edgeParam[box.edgeMap[i]] = edgeScales[i]; } std::map faceParam; int faceScales[6][2] = {{0,2},{1,2},{0,2},{1,2},{0,1},{0,1}}; for(int i = 0; i<6;i++) { - faceParam[faceLoop[i]] = &(faceScales[i]); + faceParam[box.faceMap[i]] = &(faceScales[i]); } apf::MeshIterator* it = m->begin(0); @@ -639,20 +574,20 @@ void setParameterization(gmi_model* model,apf::Mesh2* m) apf::Vector3 newParam; m->getPoint(ent,0,pt); m->getParam(ent,oldParam); - if(modelType==1 && modelTag!=sphereFaceID) + if(modelType==1 && modelTag!=sphere.faceID) { int relevantIndex = edgeParam[modelTag]; newParam[0]=pt[relevantIndex]/edgeLengths[relevantIndex]; m->setParam(ent,newParam); } - else if (modelType==2 && modelTag!=sphereFaceID && m->getDimension()>2) + else if (modelType==2 && modelTag!=sphere.faceID && m->getDimension()>2) { int* relevantIndex = faceParam[modelTag][0]; //size is 2 newParam[0] = pt[relevantIndex[0]]/edgeLengths[relevantIndex[0]]; newParam[1] = pt[relevantIndex[1]]/edgeLengths[relevantIndex[1]]; m->setParam(ent,newParam); } - else if (modelType==2 && modelTag == sphereFaceID) + else if (modelType==2 && modelTag == sphere.faceID) { double argy = (pt[1]-xyz_offset[1]); double argx = (pt[0]-xyz_offset[0]); @@ -680,17 +615,6 @@ void setParameterization(gmi_model* model,apf::Mesh2* m) m->setParam(ent,newParam); } - else if (modelType==1 && modelTag == sphereFaceID) //2D version - { - double argy = (pt[1]-xyz_offset[1]); - double argx = (pt[0]-xyz_offset[0]); - if(argx == 0 && argy ==0) - newParam[0] = 0.0; // not sure if this will happen or if this is right - else - newParam[0] = atan2(argy,argx); - - m->setParam(ent,newParam); - } } //end if } //end while @@ -698,7 +622,7 @@ void setParameterization(gmi_model* model,apf::Mesh2* m) m->acceptChanges(); } -void setParameterization2D(gmi_model* model,apf::Mesh2* m) +void setParameterization2D(gmi_model* model,apf::Mesh2* m, Enclosure box, Sphere sphere) { //Get the classification of each entity in the SimMesh apf::MeshEntity* ent; @@ -713,7 +637,7 @@ void setParameterization2D(gmi_model* model,apf::Mesh2* m) //std::cout<<"modelTag "< 139 && modelTag< 148) { - m->setModelEntity(ent,(apf::ModelEntity*)gmi_find(model,modelType,sphereFaceID)); + m->setModelEntity(ent,(apf::ModelEntity*)gmi_find(model,modelType,sphere.faceID)); } else m->setModelEntity(ent,(apf::ModelEntity*)gmi_find(model,modelType,modelTag)); @@ -731,7 +655,7 @@ void setParameterization2D(gmi_model* model,apf::Mesh2* m) double edgeLengths[2] = {boxLength,boxWidth}; for(int i=0;ibegin(0); @@ -750,14 +674,14 @@ void setParameterization2D(gmi_model* model,apf::Mesh2* m) apf::Vector3 newParam(0.0,0.0,0.0); m->getPoint(ent,0,pt); m->getParam(ent,oldParam); - if(modelType==1 && modelTag!=sphereFaceID) + if(modelType==1 && modelTag!=sphere.faceID) { int relevantIndex = edgeParam[modelTag]; newParam[0]=pt[relevantIndex]/edgeLengths[relevantIndex]; //std::cout<<"model tag "<setParam(ent,newParam); } - else if (modelType==1 && modelTag == sphereFaceID) //2D version + else if (modelType==1 && modelTag == sphere.faceID) //2D version { //std::cout<<"xyz offset "<verify(); + //add sphere - //initial adapt - initialAdapt_analytic(); + Sphere sphere = Sphere(dim); + sphere.radius = radius; - m->verify(); - return model; -} - -gmi_model* MeshAdaptPUMIDrvr::createCircleInBox(double* boxDim,double*sphereCenter, double radius) -{ - Sphere circle = Sphere(2); - circle.radius = radius; - boxLength = boxDim[0]; - boxWidth = boxDim[1]; - boxHeight = boxDim[2]; - xyz_offset[0] = sphereCenter[0]; - xyz_offset[1] = sphereCenter[1]; - xyz_offset[2] = sphereCenter[2]; - - lion_set_verbosity(1); - - //create the analytic model - gmi_model* model = gmi_make_analytic(); - - //add the sphere - - circle.makeSphere(model); + sphere.makeSphere(model); //add the box - Box box; - box.makeBox(model); - - //reparamterize circle onto box - agm_bdry b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(gmi_find(model,2,box.regionID))); - agm_use regionUse = add_adj(model, b, circle.faceID); - gmi_add_analytic_reparam(model, regionUse, reparam_Circle, 0); + geomDim = dim; + Enclosure box; + if(dim==3){ + box.makeBox3D(model); + agm_bdry b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(gmi_find(model,dim,box.regionID))); + agm_use regionUse = add_adj(model, b, sphere.faceID); + gmi_add_analytic_reparam(model, regionUse, regionFunction, 0); + + setParameterization(model,m,box,sphere); + } + else{ + box.makeBox2D(model); + agm_bdry b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(gmi_find(model,dim,box.regionID))); + agm_use regionUse = add_adj(model, b, sphere.faceID); + gmi_add_analytic_reparam(model, regionUse, reparam_Circle, 0); + setParameterization2D(model,m,box,sphere); + } - setParameterization2D(model,m); //need to modify + isAnalytic = 1; m->verify(); - //initial adapt - //initialAdapt_analytic(); - isAnalytic=1; - - m->verify(); + return; - return model; } - diff --git a/proteus/MeshAdaptPUMI/createAnalyticGeometry.h b/proteus/MeshAdaptPUMI/createAnalyticGeometry.h index a3a60b8260..4fdddd26ea 100644 --- a/proteus/MeshAdaptPUMI/createAnalyticGeometry.h +++ b/proteus/MeshAdaptPUMI/createAnalyticGeometry.h @@ -23,17 +23,4 @@ agm_bdry add_bdry(gmi_model* m, gmi_ent* e); agm_use add_adj(gmi_model* m, agm_bdry b, int tag); agm_use add_adj(gmi_model* m, agm_bdry b, int dim, int tag); -extern int edgeMap[12]; -extern int faceLoop[6]; -extern double boxLength; -extern double boxWidth; -extern double boxHeight; -extern int sphereFaceID; -extern double sphereRadius; -extern double xyz_offset[3]; - - -void makeBox(gmi_model* model); -void makeSphere(gmi_model* model); - #endif From bbb9bcad193b7cc5f9fa9f22bdc12e2f1748d76c Mon Sep 17 00:00:00 2001 From: zhang-alvin Date: Wed, 15 Jul 2020 17:49:40 -0400 Subject: [PATCH 12/16] STATE: failing to snap adapt --- proteus/MeshAdaptPUMI/AdaptHelper.py | 3 +- proteus/MeshAdaptPUMI/ErrorResidualMethod.cpp | 3 +- proteus/MeshAdaptPUMI/MeshAdaptPUMI.h | 4 +- proteus/MeshAdaptPUMI/SizeField.cpp | 2 +- proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp | 37 +++++++------------ 5 files changed, 22 insertions(+), 27 deletions(-) diff --git a/proteus/MeshAdaptPUMI/AdaptHelper.py b/proteus/MeshAdaptPUMI/AdaptHelper.py index 94d6feae87..e1d40867a2 100644 --- a/proteus/MeshAdaptPUMI/AdaptHelper.py +++ b/proteus/MeshAdaptPUMI/AdaptHelper.py @@ -662,7 +662,8 @@ def PUMI_estimateError(self): minQual = domain.AdaptManager.PUMIAdapter.getMinimumQuality() logEvent('The quality is %f ' % (minQual**(1./3.))) #adaptMeshNow=True - if(minQual**(1./3.)<0.25): + #if(minQual**(1./3.)<0.25): + if(minQual**(1./3.)<0.7): adaptMeshNow=True logEvent("Need to Adapt") diff --git a/proteus/MeshAdaptPUMI/ErrorResidualMethod.cpp b/proteus/MeshAdaptPUMI/ErrorResidualMethod.cpp index eb0992b319..8270bfdd38 100644 --- a/proteus/MeshAdaptPUMI/ErrorResidualMethod.cpp +++ b/proteus/MeshAdaptPUMI/ErrorResidualMethod.cpp @@ -944,7 +944,7 @@ void MeshAdaptPUMIDrvr::get_local_error(double &total_error) std::cerr<<"Error density maximum "<end(iter); apf::destroyField(visc); diff --git a/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h b/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h index 2180558a5f..bb197ecd04 100644 --- a/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h +++ b/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h @@ -110,13 +110,15 @@ class MeshAdaptPUMIDrvr{ bool useProteus; bool useProteusAniso; bool isAnalytic; + bool useQuality; //User Inputs std::string size_field_config; //What type of size field: interface, ERM, isotropic std::string adapt_type_config; //What type of adapt for ERM: isotropic or anisotropic - std::string logging_config; // Logging on or off + //std::string logging_config; // Logging on or off + bool logging_config; // Logging on or off //Element Residual Method void get_local_error(double& total_error); diff --git a/proteus/MeshAdaptPUMI/SizeField.cpp b/proteus/MeshAdaptPUMI/SizeField.cpp index 246d8feb53..84ee2a6a09 100644 --- a/proteus/MeshAdaptPUMI/SizeField.cpp +++ b/proteus/MeshAdaptPUMI/SizeField.cpp @@ -1444,7 +1444,7 @@ int MeshAdaptPUMIDrvr::getERMSizeField(double err_total) //apf::destroyField(curves); //apf::destroyField(hess); - if(logging_config=="on"){ + if(logging_config){ char namebuffer[20]; sprintf(namebuffer,"pumi_preadapt_aniso_%i",nAdapt); apf::writeVtkFiles(namebuffer, m); diff --git a/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp b/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp index 3ceaec3695..4ea26e962f 100644 --- a/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp +++ b/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp @@ -85,6 +85,7 @@ MeshAdaptPUMIDrvr::MeshAdaptPUMIDrvr() maxAspect = 2.0;//maxAspectRatio; hasIBM=hasInterface=hasVMS=hasERM=hasAniso=hasAnalyticSphere=useProteus=useProteusAniso=0; isAnalytic=0; + useQuality=0; } MeshAdaptPUMIDrvr::~MeshAdaptPUMIDrvr() @@ -684,7 +685,7 @@ int MeshAdaptPUMIDrvr::adaptPUMIMesh(const char* inputString) double t1 = PCU_Time(); getERMSizeField(total_error); double t2 = PCU_Time(); - if(comm_rank==0 && logging_config == "on"){ + if(comm_rank==0 && logging_config){ std::ofstream myfile; myfile.open("error_estimator_timing.txt", std::ios::app ); myfile << t2-t1<findField("proteus_sizeScale"); adapt_type_config = "anisotropic"; } -/* - if (hasTest) - testIsotropicSizeField(); -*/ -/* - if(size_field_config == "uniform"){ - //special situation where I only care about err_reg - freeField(errRho_reg); - freeField(errRel_reg); + if(useQuality){ + //size_iso=samSz::isoSize(m); + sizeFieldList.push(samSz::isoSize(m)); } -*/ + isotropicIntersect(); - if(logging_config=="on"){ + if(logging_config){ char namebuffer[50]; sprintf(namebuffer,"pumi_preadapt_%i",nAdapt); apf::writeVtkFiles(namebuffer, m); @@ -788,16 +783,8 @@ int MeshAdaptPUMIDrvr::adaptPUMIMesh(const char* inputString) in->shouldRunPostZoltan = true; in->maximumImbalance = 1.05; in->maximumIterations = numIter; -/* - if(size_field_config == "meshQuality") - { - in->shouldSnap = true; - in->shouldTransferParametric=true; - } - else - in->shouldSnap = false; -*/ in->shouldSnap=isAnalytic; + in->shouldTransferParametric=isAnalytic; //in->goodQuality = 0.16;//0.027; //double mass_before = getTotalMass(); @@ -810,7 +797,7 @@ int MeshAdaptPUMIDrvr::adaptPUMIMesh(const char* inputString) //double mass_after = getTotalMass(); //PCU_Add_Doubles(&mass_before,1); //PCU_Add_Doubles(&mass_after,1); - if(comm_rank==0 && logging_config=="on"){ + if(comm_rank==0 && logging_config){ /* std::ios::fmtflags saved(std::cout.flags()); std::cout<writeNative("DEBUG_restart.smb"); +*/ nTriggers=0; return 0; @@ -1007,6 +996,8 @@ int MeshAdaptPUMIDrvr::setAdaptProperties(std::vector sizeInputs,bo hasVMS=1; if(*it == "error_erm") hasERM=1; + if(*it == "meshQuality") + useQuality=1; } hmin=in_hmin; hmax=in_hmax; From 5bca10419c76630ae423ed0d8687f07d00eb6afe Mon Sep 17 00:00:00 2001 From: zhang-alvin Date: Tue, 28 Jul 2020 15:03:55 -0400 Subject: [PATCH 13/16] STATE: fixed parameterization at each adapt and fixed model tags leading to added mass issues in 2D --- proteus/MeshAdaptPUMI/AdaptHelper.py | 61 ++++++-- proteus/MeshAdaptPUMI/MeshAdaptPUMI.h | 3 + proteus/MeshAdaptPUMI/MeshFields.cpp | 4 + proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp | 3 + .../MeshAdaptPUMI/createAnalyticGeometry.cpp | 130 +++++++++++------- .../MeshAdaptPUMI/createAnalyticGeometry.h | 55 ++++++++ 6 files changed, 196 insertions(+), 60 deletions(-) diff --git a/proteus/MeshAdaptPUMI/AdaptHelper.py b/proteus/MeshAdaptPUMI/AdaptHelper.py index e1d40867a2..85fd644f14 100644 --- a/proteus/MeshAdaptPUMI/AdaptHelper.py +++ b/proteus/MeshAdaptPUMI/AdaptHelper.py @@ -273,7 +273,6 @@ def PUMI2Proteus(self,domain): #(cut and pasted from init, need to cleanup) self.simOutputList = [] self.auxiliaryVariables = {} - self.newAuxiliaryVariables = {} if self.simFlagsList is not None: for p, n, simFlags, model, index in zip( self.pList, @@ -318,7 +317,7 @@ def PUMI2Proteus(self,domain): av.dofsVecs = [] av.field_ids=[] av.isGaugeOwner=False - ##reinitialize auxiliaryVariables + ##reinitialize auxiliaryVariables self.auxiliaryVariables[model.name]= [av.attachModel(model,self.ar[index]) for av in n.auxiliaryVariables] else: for p,n,s,model,index in zip( @@ -335,11 +334,25 @@ def PUMI2Proteus(self,domain): for av in avList: av.attachAuxiliaryVariables(self.auxiliaryVariables) - #just added, need to initialize structures + #just added, need to initialize structures + #chrono_aux = self.auxiliaryVariables['rans2p'][0] + #chrono_aux.initialized=False + #chrono_aux.calculate_init() + #chrono_aux.model_mesh = self.modelList[1].levelModelList[-1].mesh + #chrono_aux.model_addedmass = self.modelList[2] + + #if chrono_aux.build_kdtree is True: + # chrono_aux.u = self.modelList[1].levelModelList[-1].u + # # finite element space (! linear for p, quadratic for velocity) + # chrono_aux.femSpace_velocity = chrono_aux.u[1].femSpace + # self.femSpace_pressure = chrono_aux.u[0].femSpace + # self.nodes_kdtree = spatial.cKDTree(chrono_aux.modelList[1].levelModelList[-1].mesh.nodeArray) + for model in self.modelList: logEvent("Auxiliary variable calculate_init for model %s" % (model.name,)) for av in self.auxiliaryVariables[model.name]: #gauges do not need to have aux variables recomputed as this would add additional lines to the output unnecessarily + #if hasattr(av,"isGaugeOwner") or isinstance(self.auxiliaryVariables['rans2p'][0],proteus.mbd.CouplingFSI.ProtChSystem): if hasattr(av,"isGaugeOwner"): pass else: @@ -388,6 +401,16 @@ def PUMI2Proteus(self,domain): self.modelList[self.flowIdx].levelModelList[0].coefficients.phi_s[:]=scalar[:,0] del scalar + if(self.modelList[0].name=='moveMeshElastic'): + moveModel = self.modelList[0].levelModelList[-1] + for vci in range(len(moveModel.coefficients.vectorComponents)): + moveModel.mesh.nodeVelocityArray[:,vci] = moveModel.u[vci].dof[:] + moveModel.u[vci].dof[:] = 0 + moveModel.coefficients.dt_last = modelListOld[0].levelModelList[-1].coefficients.dt_last + #import pdb; pdb.set_trace() + #if(self.modelList[2].name=='addedMass'): + # self.modelList[2].levelModelList[-1].added_mass_i = modelListOld[2].levelModelList[-1].added_mass_i + #import pdb; pdb.set_trace() logEvent("Attaching models on new mesh to each other") for m,ptmp,mOld in zip(self.modelList, self.pList, modelListOld): for lm, lu, lr, lmOld in zip(m.levelModelList, m.uList, m.rList,mOld.levelModelList): @@ -585,6 +608,17 @@ def PUMI_transferFields(self): del scalar + #mesh motion needs to be accurately transferred + #import pdb; pdb.set_trace() + if(self.modelList[0].name=='moveMeshElastic'): + vector=numpy.zeros((lm.mesh.nNodes_global,3),'d') + targetVector = self.modelList[0].levelModelList[-1].mesh.nodeVelocityArray + for vci in range(targetVector.shape[1]): + vector[:,vci] = self.modelList[0].levelModelList[0].mesh.nodeVelocityArray[:,vci] + domain.AdaptManager.PUMIAdapter.transferFieldToPUMI( + self.modelList[0].levelModelList[0].coefficients.vectorName.encode('utf-8'), vector) + del vector + #Get Physical Parameters #Can we do this in a problem-independent way? @@ -618,6 +652,18 @@ def PUMI_estimateError(self): domain.AdaptManager.PUMIAdapter.adaptMesh() and self.so.useOneMesh): #and #self.nSolveSteps%domain.AdaptManager.PUMIAdapter.numAdaptSteps()==0): + + #this needs to update before mesh coordinates are updated for reparameterization + if (self.auxiliaryVariables['rans2p'][0].subcomponents[0].__class__.__name__== 'ProtChBody'): + sphereCoords = numpy.asarray(self.auxiliaryVariables['rans2p'][0].subcomponents[0].position) + logEvent("The sphere coordinates are %f %f %f" % (sphereCoords[0],sphereCoords[1],sphereCoords[2])) + domain.AdaptManager.PUMIAdapter.updateSphereCoordinates(sphereCoords) + logEvent("Updated the sphere coordinates %f %f %f" % (sphereCoords[0],sphereCoords[1],sphereCoords[2])) + else: + sys.exit("Haven't been implemented code yet to cover this behavior.") + + + if (b"isotropicProteus" in self.domain.AdaptManager.sizeInputs): domain.AdaptManager.PUMIAdapter.transferFieldToPUMI("proteus_size", self.modelList[0].levelModelList[0].mesh.size_field) @@ -667,16 +713,9 @@ def PUMI_estimateError(self): adaptMeshNow=True logEvent("Need to Adapt") - if (self.auxiliaryVariables['rans2p'][0].subcomponents[0].__class__.__name__== 'ProtChBody'): - sphereCoords = numpy.asarray(self.auxiliaryVariables['rans2p'][0].subcomponents[0].position) - domain.AdaptManager.PUMIAdapter.updateSphereCoordinates(sphereCoords) - logEvent("Updated the sphere coordinates %f %f %f" % (sphereCoords[0],sphereCoords[1],sphereCoords[2])) - else: - sys.exit("Haven't been implemented code yet to cover this behavior.") - if(b'pseudo' in self.domain.AdaptManager.sizeInputs): adaptMeshNow=True - + #if not adapting need to return data structures to original form which was modified by PUMI_transferFields() if(adaptMeshNow == False): for m in self.modelList: diff --git a/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h b/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h index bb197ecd04..033de9fdff 100644 --- a/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h +++ b/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h @@ -4,6 +4,7 @@ #include #include #include "PyEmbeddedFunctions.h" +#include "createAnalyticGeometry.h" /** \file MeshAdaptPUMI.h @@ -74,6 +75,8 @@ class MeshAdaptPUMIDrvr{ //analytic geometry void createAnalyticGeometry(int dim, double* boxDim, double*sphereCenter,double radius); + Enclosure modelBox; + Sphere modelSphere; void updateSphereCoordinates(double*sphereCenter); void initialAdapt_analytic(); diff --git a/proteus/MeshAdaptPUMI/MeshFields.cpp b/proteus/MeshAdaptPUMI/MeshFields.cpp index 2d18cde5e2..c7fa1610d6 100644 --- a/proteus/MeshAdaptPUMI/MeshFields.cpp +++ b/proteus/MeshAdaptPUMI/MeshFields.cpp @@ -80,6 +80,10 @@ int MeshAdaptPUMIDrvr::transferFieldToPUMI(const char* name, double const* inArr apf::setComponents(f, v, 0, &tmp[0]); } m->end(it); + if (!strcmp(name, "coordinates") && isAnalytic) { + Reparam::reparameterizeEntities(m->getModel(),m,modelBox,modelSphere); + } + return 0; } diff --git a/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp b/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp index 4ea26e962f..ed72226988 100644 --- a/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp +++ b/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp @@ -783,6 +783,8 @@ int MeshAdaptPUMIDrvr::adaptPUMIMesh(const char* inputString) in->shouldRunPostZoltan = true; in->maximumImbalance = 1.05; in->maximumIterations = numIter; + //in->shouldSnap=0; + //in->shouldTransferParametric=0; in->shouldSnap=isAnalytic; in->shouldTransferParametric=isAnalytic; //in->goodQuality = 0.16;//0.027; @@ -836,6 +838,7 @@ int MeshAdaptPUMIDrvr::adaptPUMIMesh(const char* inputString) */ nTriggers=0; + //std::exit(1); return 0; } diff --git a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp index 6dffe6476a..b97833dc86 100644 --- a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp +++ b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp @@ -20,19 +20,18 @@ agm_use add_adj(gmi_model* m, agm_bdry b, int tag) //anonymous namespace to make all functions local to this file scope -namespace -{ - int edgeMap[12] = {50,48,46,52,11,16,20,6,73,72,71,74}; - int faceLoop[6] = {80,78,76,82,42,24}; - - double sphereRadius; - double xyz_offset[3]; +double xyz_offset[3]; +double sphereRadius; +double boxLength; +double boxWidth; +double boxHeight; +int geomDim; +Enclosure modelBox; +Sphere modelSphere; - double boxLength; - double boxWidth; - double boxHeight; - int geomDim; +namespace +{ void vert0(double const p[2], double x[3], void*) { (void)p; @@ -298,24 +297,6 @@ namespace typedef void (*EntityMapArray) (double const p[2], double x[3], void*); typedef void (*ParametricFunctionArray) (double const from[2], double to[2], void*); -struct Enclosure{ - - std::vector vertexMap; - double vertRan[1][2]={{0.0,0.0}}; - int vertPer=0; - - std::vector edgeMap; - int edgePer = 0; - double edgeRan[1][2] = {{0.0,1.0}}; - - std::vector faceMap; - int faPer[2] = {0, 0}; - double faRan[2][2] = {{0,1},{0,1}}; - - int regionID = 92; // fixed ID - void makeBox2D(gmi_model* model); - void makeBox3D(gmi_model* model); -}; void Enclosure::makeBox2D(gmi_model* model) { @@ -333,7 +314,7 @@ void Enclosure::makeBox2D(gmi_model* model) for(auto i=0; iacceptChanges(); } +void Reparam::reparameterizeEntities(gmi_model*model,apf::Mesh2*m,Enclosure box, Sphere sphere){ + std::map edgeParam; + const int numEdges = 4; + int edgeScales[numEdges] = {0,1,0,1}; + double edgeLengths[2] = {boxLength,boxWidth}; + for(int i=0;ibegin(0); + apf::MeshEntity* ent; + while( (ent = m->iterate(it)) ) + { + apf::ModelEntity* g_ent = m->toModel(ent); + + apf::MeshEntity* ev[2]; + m->getDownward(ent,0,ev); + int modelTag = m->getModelTag(g_ent); + int modelType = m->getModelType(g_ent); + if(modelType<3 && modelType!=0) + { + apf::Vector3 pt; + apf::Vector3 oldParam; + apf::Vector3 newParam(0.0,0.0,0.0); + m->getPoint(ent,0,pt); + m->getParam(ent,oldParam); + if(modelType==1 && modelTag!=sphere.faceID) + { + int relevantIndex = edgeParam[modelTag]; + newParam[0]=pt[relevantIndex]/edgeLengths[relevantIndex]; + //std::cout<<"model tag "<setParam(ent,newParam); + } + else if (modelType==1 && modelTag == sphere.faceID) //2D version + { + //std::cout<<"xyz offset "<setParam(ent,newParam); + } + else if (modelType==2) + { + newParam[0] = pt[0]/boxLength; + newParam[1] = pt[1]/boxWidth; + m->setParam(ent,newParam); + } + + } //end if + } //end while + m->end(it); + m->acceptChanges(); + +} + void setParameterization2D(gmi_model* model,apf::Mesh2* m, Enclosure box, Sphere sphere) { //Get the classification of each entity in the SimMesh @@ -649,6 +671,7 @@ void setParameterization2D(gmi_model* model,apf::Mesh2* m, Enclosure box, Sphere m->acceptChanges(); //Need to set the parametric coordinates of each of the boundary vertices +/* std::map edgeParam; const int numEdges = 4; int edgeScales[numEdges] = {0,1,0,1}; @@ -704,6 +727,8 @@ void setParameterization2D(gmi_model* model,apf::Mesh2* m, Enclosure box, Sphere } //end while m->end(it); m->acceptChanges(); +*/ + Reparam::reparameterizeEntities(model,m,box,sphere); } void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ @@ -784,6 +809,11 @@ void MeshAdaptPUMIDrvr::updateSphereCoordinates(double*sphereCenter) xyz_offset[0] = sphereCenter[0]; xyz_offset[1] = sphereCenter[1]; xyz_offset[2] = sphereCenter[2]; + + char buffer[100]; + sprintf(buffer,"Checking coordinates at update %f %f %f",xyz_offset[0],xyz_offset[1],xyz_offset[2]); + logEvent(buffer,3); + } void MeshAdaptPUMIDrvr::createAnalyticGeometry(int dim, double* boxDim,double*sphereCenter, double radius) @@ -822,6 +852,8 @@ void MeshAdaptPUMIDrvr::createAnalyticGeometry(int dim, double* boxDim,double*sp setParameterization2D(model,m,box,sphere); } + modelBox = box; + modelSphere = sphere; isAnalytic = 1; m->verify(); diff --git a/proteus/MeshAdaptPUMI/createAnalyticGeometry.h b/proteus/MeshAdaptPUMI/createAnalyticGeometry.h index 4fdddd26ea..3f8f276643 100644 --- a/proteus/MeshAdaptPUMI/createAnalyticGeometry.h +++ b/proteus/MeshAdaptPUMI/createAnalyticGeometry.h @@ -9,6 +9,7 @@ #include #include #include +#include #include #include @@ -23,4 +24,58 @@ agm_bdry add_bdry(gmi_model* m, gmi_ent* e); agm_use add_adj(gmi_model* m, agm_bdry b, int tag); agm_use add_adj(gmi_model* m, agm_bdry b, int dim, int tag); +extern double boxLength; +extern double boxWidth; +extern double boxHeight; +extern double sphereRadius; +extern double xyz_offset[3]; + +struct Enclosure{ + + std::vector vertexMap; + double vertRan[1][2]={{0.0,0.0}}; + int vertPer=0; + + std::vector edgeMap; + int edgePer = 0; + double edgeRan[1][2] = {{0.0,1.0}}; + + std::vector faceMap; + int faPer[2] = {0, 0}; + double faRan[2][2] = {{0,1},{0,1}}; + + int regionID = 92; // fixed ID + void makeBox2D(gmi_model* model); + void makeBox3D(gmi_model* model); +}; + +class Sphere{ + public: + const double pi = apf::pi; + int faceID = 6;//123; + double radius; + double offset[3]; + int dim; + int faPer[2] = {1, 0}; + double faRan[2][2] = {{0,6.28318530718},{0.0,apf::pi}}; + + Sphere(int x){ + dim = x; + if(dim==2){ + faRan[1][1] = 0.0; + } + } + Sphere(){} + void makeSphere(gmi_model* model); +}; + +namespace Reparam{ + void reparameterizeEntities(gmi_model*model,apf::Mesh2*m,Enclosure box, Sphere sphere); +} + +//namespace MeshAdaptPUMIDrvr{ +//extern Enclosure modelBox; +//extern Sphere modelSphere; +//} + #endif From 697ce29db54e3de517e1c1b4a2f41d7e55514db6 Mon Sep 17 00:00:00 2001 From: zhang-alvin Date: Fri, 31 Jul 2020 17:53:30 -0400 Subject: [PATCH 14/16] STATE: working box only in 3D. had to revert to previous setup because parameterization is tricky --- proteus/MeshAdaptPUMI/AdaptHelper.py | 15 +- .../MeshAdaptPUMI/createAnalyticGeometry.cpp | 236 ++++++++++++------ .../MeshAdaptPUMI/createAnalyticGeometry.h | 13 +- 3 files changed, 175 insertions(+), 89 deletions(-) diff --git a/proteus/MeshAdaptPUMI/AdaptHelper.py b/proteus/MeshAdaptPUMI/AdaptHelper.py index 85fd644f14..6c4a0a29c4 100644 --- a/proteus/MeshAdaptPUMI/AdaptHelper.py +++ b/proteus/MeshAdaptPUMI/AdaptHelper.py @@ -654,13 +654,14 @@ def PUMI_estimateError(self): #self.nSolveSteps%domain.AdaptManager.PUMIAdapter.numAdaptSteps()==0): #this needs to update before mesh coordinates are updated for reparameterization - if (self.auxiliaryVariables['rans2p'][0].subcomponents[0].__class__.__name__== 'ProtChBody'): - sphereCoords = numpy.asarray(self.auxiliaryVariables['rans2p'][0].subcomponents[0].position) - logEvent("The sphere coordinates are %f %f %f" % (sphereCoords[0],sphereCoords[1],sphereCoords[2])) - domain.AdaptManager.PUMIAdapter.updateSphereCoordinates(sphereCoords) - logEvent("Updated the sphere coordinates %f %f %f" % (sphereCoords[0],sphereCoords[1],sphereCoords[2])) - else: - sys.exit("Haven't been implemented code yet to cover this behavior.") + if self.auxiliaryVariables['rans2p'] != []: + if(self.auxiliaryVariables['rans2p'][0].subcomponents[0].__class__.__name__== 'ProtChBody'): + sphereCoords = numpy.asarray(self.auxiliaryVariables['rans2p'][0].subcomponents[0].position) + logEvent("The sphere coordinates are %f %f %f" % (sphereCoords[0],sphereCoords[1],sphereCoords[2])) + domain.AdaptManager.PUMIAdapter.updateSphereCoordinates(sphereCoords) + logEvent("Updated the sphere coordinates %f %f %f" % (sphereCoords[0],sphereCoords[1],sphereCoords[2])) + else: + sys.exit("Haven't been implemented code yet to cover this behavior.") diff --git a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp index b97833dc86..3bb27ec561 100644 --- a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp +++ b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp @@ -257,6 +257,28 @@ namespace to[1] = from[0]; } + void reparamREdge_0(double const from[2], double to[2], void*) + { + to[0] = from[0]; + to[1] = 0.0; + } + + void reparamREdge_1(double const from[2], double to[2], void*) + { + to[0] = 0.0; + to[1] = from[0]; + } + void reparamREdge_2(double const from[2], double to[2], void*) + { + to[0] = from[0]; + to[1] = 1.0; + } + + void reparamREdge_3(double const from[2], double to[2], void*) + { + to[0] = 1.0; + to[1] = from[0]; + } void regionFunction(double const p[2], double x[3], void*) { @@ -403,6 +425,11 @@ void Enclosure::makeBox3D(gmi_model* model) std::vector> indexPairs = {{0,1}, {1,2}, {2,3}, {3,0}, {4,5},{5,6},{6,7},{7,4},{0,4},{1,5},{2,6},{3,7} }; + + //std::vector> indexPairs = {{0,1}, {1,2}, {3,2}, {0,3}, + // {4,5},{5,6},{7,6},{4,7},{0,4},{1,5},{2,6},{3,7} }; + +/* for(int i=0;isetModelEntity(ent,(apf::ModelEntity*)gmi_find(model,modelType,sphere.faceID)); } - else + else{ m->setModelEntity(ent,(apf::ModelEntity*)gmi_find(model,modelType,modelTag)); + } } m->end(it); } m->setModel(model); m->acceptChanges(); - + + Reparam::reparameterizeEntities(model,m,box,sphere); + +} + +void Reparam::reparameterize3D(gmi_model*model,apf::Mesh2*m,Enclosure box, Sphere sphere){ //Need to set the parametric coordinates of each of the boundary vertices std::map edgeParam; @@ -521,6 +630,7 @@ void setParameterization(gmi_model* model,apf::Mesh2* m, Enclosure box, Sphere s } apf::MeshIterator* it = m->begin(0); + apf::MeshEntity* ent; while( (ent = m->iterate(it)) ) { apf::ModelEntity* g_ent = m->toModel(ent); @@ -532,17 +642,16 @@ void setParameterization(gmi_model* model,apf::Mesh2* m, Enclosure box, Sphere s if(modelType<3 && modelType!=0) { apf::Vector3 pt; - apf::Vector3 oldParam; apf::Vector3 newParam; m->getPoint(ent,0,pt); - m->getParam(ent,oldParam); - if(modelType==1 && modelTag!=sphere.faceID) + if(modelType==1) { int relevantIndex = edgeParam[modelTag]; newParam[0]=pt[relevantIndex]/edgeLengths[relevantIndex]; m->setParam(ent,newParam); + std::cout<<"modelType? "<getDimension()>2) + else if (modelType==2 && modelTag!=sphere.faceID) { int* relevantIndex = faceParam[modelTag][0]; //size is 2 newParam[0] = pt[relevantIndex[0]]/edgeLengths[relevantIndex[0]]; @@ -582,9 +691,17 @@ void setParameterization(gmi_model* model,apf::Mesh2* m, Enclosure box, Sphere s } //end while m->end(it); m->acceptChanges(); + } void Reparam::reparameterizeEntities(gmi_model*model,apf::Mesh2*m,Enclosure box, Sphere sphere){ + if(m->getDimension()==2) + reparameterize2D(model,m,box, sphere); + else + reparameterize3D(model,m,box, sphere); +} + +void Reparam::reparameterize2D(gmi_model*model,apf::Mesh2*m,Enclosure box, Sphere sphere){ std::map edgeParam; const int numEdges = 4; int edgeScales[numEdges] = {0,1,0,1}; @@ -671,69 +788,13 @@ void setParameterization2D(gmi_model* model,apf::Mesh2* m, Enclosure box, Sphere m->acceptChanges(); //Need to set the parametric coordinates of each of the boundary vertices -/* - std::map edgeParam; - const int numEdges = 4; - int edgeScales[numEdges] = {0,1,0,1}; - double edgeLengths[2] = {boxLength,boxWidth}; - for(int i=0;ibegin(0); - while( (ent = m->iterate(it)) ) - { - apf::ModelEntity* g_ent = m->toModel(ent); - - apf::MeshEntity* ev[2]; - m->getDownward(ent,0,ev); - int modelTag = m->getModelTag(g_ent); - int modelType = m->getModelType(g_ent); - if(modelType<3 && modelType!=0) - { - apf::Vector3 pt; - apf::Vector3 oldParam; - apf::Vector3 newParam(0.0,0.0,0.0); - m->getPoint(ent,0,pt); - m->getParam(ent,oldParam); - if(modelType==1 && modelTag!=sphere.faceID) - { - int relevantIndex = edgeParam[modelTag]; - newParam[0]=pt[relevantIndex]/edgeLengths[relevantIndex]; - //std::cout<<"model tag "<setParam(ent,newParam); - } - else if (modelType==1 && modelTag == sphere.faceID) //2D version - { - //std::cout<<"xyz offset "<setParam(ent,newParam); - } - else if (modelType==2) - { - newParam[0] = pt[0]/boxLength; - newParam[1] = pt[1]/boxWidth; - m->setParam(ent,newParam); - } - - } //end if - } //end while - m->end(it); - m->acceptChanges(); -*/ - Reparam::reparameterizeEntities(model,m,box,sphere); + Reparam::reparameterizeEntities(model,m,box,sphere); } void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ //at this point, hmin and hmax haven't been set yet //apf::Field* size_initial = apf::createLagrangeField(m,"size_initial",apf::SCALAR,1); + lion_set_verbosity(1); size_iso = apf::createLagrangeField(m,"proteus_size",apf::SCALAR,1); apf::MeshIterator* it = m->begin(0); apf::MeshEntity* ent; @@ -757,6 +818,28 @@ void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ //gradeMesh(1.5); + apf::Field* modelTagge = apf::createLagrangeField(m,"modelTags",apf::SCALAR,1); + it = m->begin(0); + while( (ent = m->iterate(it)) ) + { + apf::ModelEntity* g_ent = m->toModel(ent); + int modelTag = m->getModelTag(g_ent); + int modelType = m->getModelType(g_ent); + apf::setScalar(modelTagge,ent,0,modelTag); + } + m->end(it); + + apf::Field* modelTyper = apf::createLagrangeField(m,"modelType",apf::SCALAR,1); + it = m->begin(0); + while( (ent = m->iterate(it)) ) + { + apf::ModelEntity* g_ent = m->toModel(ent); + int modelTag = m->getModelTag(g_ent); + int modelType = m->getModelType(g_ent); + apf::setScalar(modelTyper,ent,0,modelType); + } + m->end(it); + apf::writeVtkFiles("initialProteus",m); ma::Input* in = ma::configure(m,size_iso); in->maximumIterations = 10; @@ -764,7 +847,7 @@ void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ in->shouldTransferParametric = true; in->shouldFixShape = true; in->debugFolder="./debug_fine"; - ma::adaptVerbose(in,false); + ma::adaptVerbose(in,true); m->verify(); apf::writeVtkFiles("middleProteus",m); @@ -783,7 +866,7 @@ void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ else apf::setScalar(size_iso,ent,0,hmax); */ - apf::setScalar(size_iso,ent,0,0.01); + apf::setScalar(size_iso,ent,0,1.0); } m->end(it); @@ -838,10 +921,9 @@ void MeshAdaptPUMIDrvr::createAnalyticGeometry(int dim, double* boxDim,double*sp Enclosure box; if(dim==3){ box.makeBox3D(model); - agm_bdry b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(gmi_find(model,dim,box.regionID))); - agm_use regionUse = add_adj(model, b, sphere.faceID); - gmi_add_analytic_reparam(model, regionUse, regionFunction, 0); - + //agm_bdry b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(gmi_find(model,dim,box.regionID))); + //agm_use regionUse = add_adj(model, b, sphere.faceID); + //gmi_add_analytic_reparam(model, regionUse, regionFunction, 0); setParameterization(model,m,box,sphere); } else{ @@ -857,6 +939,10 @@ void MeshAdaptPUMIDrvr::createAnalyticGeometry(int dim, double* boxDim,double*sp isAnalytic = 1; m->verify(); + //initialAdapt_analytic(); + //initialAdapt_analytic(1.0); + //initialAdapt_analytic(0.1); + return; } diff --git a/proteus/MeshAdaptPUMI/createAnalyticGeometry.h b/proteus/MeshAdaptPUMI/createAnalyticGeometry.h index 3f8f276643..cc861fce6a 100644 --- a/proteus/MeshAdaptPUMI/createAnalyticGeometry.h +++ b/proteus/MeshAdaptPUMI/createAnalyticGeometry.h @@ -51,8 +51,7 @@ struct Enclosure{ class Sphere{ public: - const double pi = apf::pi; - int faceID = 6;//123; + int faceID; //values are dictated by spatial tools double radius; double offset[3]; int dim; @@ -63,7 +62,10 @@ class Sphere{ dim = x; if(dim==2){ faRan[1][1] = 0.0; + faceID=6; } + if(dim==3) + faceID=9; } Sphere(){} void makeSphere(gmi_model* model); @@ -71,11 +73,8 @@ class Sphere{ namespace Reparam{ void reparameterizeEntities(gmi_model*model,apf::Mesh2*m,Enclosure box, Sphere sphere); + void reparameterize2D(gmi_model*model,apf::Mesh2*m,Enclosure box, Sphere sphere); + void reparameterize3D(gmi_model*model,apf::Mesh2*m,Enclosure box, Sphere sphere); } -//namespace MeshAdaptPUMIDrvr{ -//extern Enclosure modelBox; -//extern Sphere modelSphere; -//} - #endif From 601dee541c670adad6615fd055cd561135664430 Mon Sep 17 00:00:00 2001 From: zhang-alvin Date: Thu, 6 Aug 2020 20:02:44 -0400 Subject: [PATCH 15/16] ENH: added in memory splitting/partitioning --- .../MeshAdaptPUMI/createAnalyticGeometry.cpp | 74 +++++++++++++++++-- .../MeshAdaptPUMI/createAnalyticGeometry.h | 36 +++++++++ 2 files changed, 104 insertions(+), 6 deletions(-) diff --git a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp index 3bb27ec561..6e53a80875 100644 --- a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp +++ b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp @@ -649,7 +649,7 @@ void Reparam::reparameterize3D(gmi_model*model,apf::Mesh2*m,Enclosure box, Spher int relevantIndex = edgeParam[modelTag]; newParam[0]=pt[relevantIndex]/edgeLengths[relevantIndex]; m->setParam(ent,newParam); - std::cout<<"modelType? "<end(it); @@ -866,7 +866,7 @@ void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ else apf::setScalar(size_iso,ent,0,hmax); */ - apf::setScalar(size_iso,ent,0,1.0); + apf::setScalar(size_iso,ent,0,0.25); } m->end(it); @@ -899,6 +899,7 @@ void MeshAdaptPUMIDrvr::updateSphereCoordinates(double*sphereCenter) } +int Splitter::partitionFactor; void MeshAdaptPUMIDrvr::createAnalyticGeometry(int dim, double* boxDim,double*sphereCenter, double radius) { boxLength = boxDim[0]; @@ -921,9 +922,9 @@ void MeshAdaptPUMIDrvr::createAnalyticGeometry(int dim, double* boxDim,double*sp Enclosure box; if(dim==3){ box.makeBox3D(model); - //agm_bdry b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(gmi_find(model,dim,box.regionID))); - //agm_use regionUse = add_adj(model, b, sphere.faceID); - //gmi_add_analytic_reparam(model, regionUse, regionFunction, 0); + agm_bdry b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(gmi_find(model,dim,box.regionID))); + agm_use regionUse = add_adj(model, b, sphere.faceID); + gmi_add_analytic_reparam(model, regionUse, regionFunction, 0); setParameterization(model,m,box,sphere); } else{ @@ -939,6 +940,20 @@ void MeshAdaptPUMIDrvr::createAnalyticGeometry(int dim, double* boxDim,double*sp isAnalytic = 1; m->verify(); + //m->writeNative("analyticMesh.smb"); + apf::Migration* plan = 0; + gmi_register_mesh(); + Splitter::partitionFactor = PCU_Comm_Peers(); + Splitter::switchToOriginals(); + bool isOriginal = ((PCU_Comm_Self() % Splitter::partitionFactor) == 0); + if (isOriginal) { + plan = Splitter::getPlan(m); + } + Splitter::switchToAll(); + m = apf::repeatMdsMesh(m, model, plan, Splitter::partitionFactor); + Parma_PrintPtnStats(m, ""); + //m->writeNative("analyticMesh.smb"); + //initialAdapt_analytic(); //initialAdapt_analytic(1.0); //initialAdapt_analytic(0.1); @@ -946,3 +961,50 @@ void MeshAdaptPUMIDrvr::createAnalyticGeometry(int dim, double* boxDim,double*sp return; } + + +void Splitter::freeMesh(apf::Mesh* m) +{ + m->destroyNative(); + apf::destroyMesh(m); +} + +apf::Migration* Splitter::getPlan(apf::Mesh* m) +{ +/* + apf::Splitter* splitter = apf::makeZoltanSplitter( + m, apf::GRAPH, apf::PARTITION, false); + apf::MeshTag* weights = Parma_WeighByMemory(m); + apf::Migration* plan = splitter->split(weights, 1.05, partitionFactor); +*/ + + apf::Splitter* splitter = Parma_MakeRibSplitter(m); + apf::MeshTag* weights = Parma_WeighByMemory(m); + apf::Migration* plan = splitter->split(weights, 1.10, partitionFactor); + + + apf::removeTagFromDimension(m, weights, m->getDimension()); + m->destroyTag(weights); + delete splitter; + return plan; +} + +void Splitter::switchToOriginals() +{ + int self = PCU_Comm_Self(); + int groupRank = self / Splitter::partitionFactor; + int group = self % Splitter::partitionFactor; + MPI_Comm groupComm; + MPI_Comm_split(MPI_COMM_WORLD, group, groupRank, &groupComm); + PCU_Switch_Comm(groupComm); +} + +void Splitter::switchToAll() +{ + MPI_Comm prevComm = PCU_Get_Comm(); + PCU_Switch_Comm(MPI_COMM_WORLD); + MPI_Comm_free(&prevComm); + PCU_Barrier(); +} + + diff --git a/proteus/MeshAdaptPUMI/createAnalyticGeometry.h b/proteus/MeshAdaptPUMI/createAnalyticGeometry.h index cc861fce6a..7b0241c8c9 100644 --- a/proteus/MeshAdaptPUMI/createAnalyticGeometry.h +++ b/proteus/MeshAdaptPUMI/createAnalyticGeometry.h @@ -77,4 +77,40 @@ namespace Reparam{ void reparameterize3D(gmi_model*model,apf::Mesh2*m,Enclosure box, Sphere sphere); } +#include +#include +#include +#include + +namespace Splitter{ + + extern int partitionFactor; + void freeMesh(apf::Mesh* m); + apf::Migration* getPlan(apf::Mesh* m); + void switchToOriginals(); + void switchToAll(); + +} +/* +int main(int argc, char** argv) +{ + gmi_register_mesh(); + getConfig(argc,argv); + bool isOriginal = ((PCU_Comm_Self() % partitionFactor) == 0); + gmi_model* g = 0; + g = gmi_load(modelFile); + apf::Mesh2* m = 0; + apf::Migration* plan = 0; + switchToOriginals(); + if (isOriginal) { + m = apf::loadMdsMesh(g, meshFile); + plan = getPlan(m); + } + switchToAll(); + m = apf::repeatMdsMesh(m, g, plan, partitionFactor); + Parma_PrintPtnStats(m, ""); + m->writeNative(outFile); + freeMesh(m); +*/ + #endif From 6f76661ffe4abfb503c3171aae8f306929e3ce53 Mon Sep 17 00:00:00 2001 From: zhang-alvin Date: Wed, 19 Aug 2020 11:35:34 -0400 Subject: [PATCH 16/16] STATE: piercing cylinder, but not functional 3D adapt --- proteus/MeshAdaptPUMI/MeshAdaptPUMI.h | 4 + proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp | 33 +- .../MeshAdaptPUMI/createAnalyticGeometry.cpp | 406 ++++++++++++++++-- .../MeshAdaptPUMI/createAnalyticGeometry.h | 18 + 4 files changed, 426 insertions(+), 35 deletions(-) diff --git a/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h b/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h index 033de9fdff..c2149d332a 100644 --- a/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h +++ b/proteus/MeshAdaptPUMI/MeshAdaptPUMI.h @@ -75,8 +75,12 @@ class MeshAdaptPUMIDrvr{ //analytic geometry void createAnalyticGeometry(int dim, double* boxDim, double*sphereCenter,double radius); + void createAnalyticGeometryCylinder(int dim, double* boxDim, double*sphereCenter,double radius); Enclosure modelBox; Sphere modelSphere; + Sphere modelCircle1; + Sphere modelCircle2; + PiercingCylinder modelPiercingCylinder; void updateSphereCoordinates(double*sphereCenter); void initialAdapt_analytic(); diff --git a/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp b/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp index ed72226988..8fba9a2081 100644 --- a/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp +++ b/proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp @@ -167,8 +167,13 @@ int MeshAdaptPUMIDrvr::loadMeshForAnalytic(const char* meshFile,double* boxDim,d m->verify(); //create analytic geometry - createAnalyticGeometry(m->getDimension(),boxDim,sphereCenter,radius); +// createAnalyticGeometry(m->getDimension(),boxDim,sphereCenter,radius); + createAnalyticGeometryCylinder(m->getDimension(),boxDim,sphereCenter,radius); + //m->writeNative("analyticMesh.smb"); + //initialAdapt_analytic(); + //std::exit(1); m->verify(); + return 0; } @@ -718,7 +723,26 @@ int MeshAdaptPUMIDrvr::adaptPUMIMesh(const char* inputString) } if(useQuality){ //size_iso=samSz::isoSize(m); - sizeFieldList.push(samSz::isoSize(m)); + apf::MeshIterator* it = m->begin(0); + apf::MeshEntity* ent; + + //sizeFieldList.push(samSz::isoSize(m)); + apf::Field * test = m->findField("proteus_init"); + freeField(test); + apf::Field* size_init = apf::createLagrangeField(m,"proteus_init",apf::SCALAR,1); + while( (ent = m->iterate(it)) ) + { + apf::ModelEntity* g_ent = m->toModel(ent); + int modelTag = m->getModelTag(g_ent); + int modelType = m->getModelType(g_ent); + if(modelTag == modelPiercingCylinder.faceID && modelType==2 || modelType==1 && (modelTag == modelCircle1.faceID || modelTag == modelCircle2.faceID)) + apf::setScalar(size_init,ent,0,0.05); + //apf::setScalar(size_init,ent,0,0.2); + else + apf::setScalar(size_init,ent,0,0.2); + } + m->end(it); + sizeFieldList.push(size_init); } isotropicIntersect(); @@ -792,7 +816,8 @@ int MeshAdaptPUMIDrvr::adaptPUMIMesh(const char* inputString) double t1 = PCU_Time(); //ma::adapt(in); - ma::adaptVerbose(in); + in->debugFolder="./debug_fine"; + ma::adaptVerbose(in,true); double t2 = PCU_Time(); m->verify(); @@ -1017,3 +1042,5 @@ int MeshAdaptPUMIDrvr::setAdaptProperties(std::vector sizeInputs,bo return 0; } + + diff --git a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp index 6e53a80875..a5c03b4bf1 100644 --- a/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp +++ b/proteus/MeshAdaptPUMI/createAnalyticGeometry.cpp @@ -28,7 +28,9 @@ double boxHeight; int geomDim; Enclosure modelBox; Sphere modelSphere; - +Sphere modelCircle1; +Sphere modelCircle2; +PiercingCylinder modelPiercingCylinder; namespace { @@ -295,6 +297,48 @@ namespace to[0] = x/boxLength; to[1] = y/boxWidth; } + void reparam_Circle0(double const from[2], double to[2], void*){ + + //given theta, need to get y and x in parameterized form + double x = sphereRadius*cos(from[0])+xyz_offset[0]; + double y = sphereRadius*sin(from[0])+xyz_offset[1]; + to[0] = x/boxLength; + to[1] = y/boxWidth; + to[2] = 0.0; + } + void reparam_Circle1(double const from[2], double to[2], void*){ + + //given theta, need to get y and x in parameterized form + double x = sphereRadius*cos(from[0])+xyz_offset[0]; + double y = sphereRadius*sin(from[0])+xyz_offset[1]; + to[0] = x/boxLength; + to[1] = y/boxWidth; + to[2] = boxHeight; + } + + void reparam_CircleCylinder0(double const from[2], double to[2], void*){ + + //given theta, need to get y and x in parameterized form + //double x = sphereRadius*cos(from[0])+xyz_offset[0]; + //double y = sphereRadius*sin(from[0])+xyz_offset[1]; + //to[0] = x/boxLength; + //to[1] = y/boxWidth; + to[0] = from[0]; + to[1] = 0.0;//from[1]; + //to[2] = 0; + } + void reparam_CircleCylinder1(double const from[2], double to[2], void*){ + + //given theta, need to get y and x in parameterized form + //double x = sphereRadius*cos(from[0])+xyz_offset[0]; + //double y = sphereRadius*sin(from[0])+xyz_offset[1]; + //to[0] = x/boxLength; + //to[1] = y/boxWidth; + to[0] = from[0]; + to[1] = 1.0;//from[1]; + //to[2] = 1.0;//boxHeight; + } + //need to set these functions separately from a member because the function signatures are otherwise modified. //Likewise the internals need to use static variables to still connect them with other objects @@ -312,6 +356,74 @@ namespace x[2] = xyz_offset[2]+sphereRadius*cos(p[1]); } } + void circleFace0(double const p[2], double x[3], void*) + { + x[0] = xyz_offset[0]+sphereRadius*cos(p[0]); + x[1] = xyz_offset[1]+sphereRadius*sin(p[0]); + x[2] = 0.0; + } + + void circleFace1(double const p[2], double x[3], void*) + { + x[0] = xyz_offset[0]+sphereRadius*cos(p[0]); + x[1] = xyz_offset[1]+sphereRadius*sin(p[0]); + x[2] = boxHeight; + } + + void cylinderFace(double const p[2], double x[3], void*) + { + x[0] = xyz_offset[0]+sphereRadius*cos(p[0]); + //x[1] = xyz_offset[1]+sphereRadius*sin(p[0]) * sin(p[1]); + x[1] = xyz_offset[1]+sphereRadius*sin(p[0]); + x[2] = boxHeight*p[1]; //need to double check this + } + + +} + +void checkEntities(apf::Mesh* m){ + apf::MeshEntity* ent; + apf::MeshIterator* it; +/* + for(int i =0;i<4;i++) + { + apf::MeshIterator* it = m->begin(i); + while( (ent = m->iterate(it))) + { + apf::ModelEntity* g_ent = m->toModel(ent); + int modelTag = m->getModelTag(g_ent); + int modelType = m->getModelType(g_ent); + std::cout<<"model Tag, type "<end(it); + } +*/ + if(m->findField("modelTags")) + apf::destroyField(m->findField("modelTags")); + if(m->findField("modelType")) + apf::destroyField(m->findField("modelType")); + + apf::Field* modelTagge = apf::createLagrangeField(m,"modelTags",apf::SCALAR,1); + it = m->begin(0); + while( (ent = m->iterate(it)) ) + { + apf::ModelEntity* g_ent = m->toModel(ent); + int modelTag = m->getModelTag(g_ent); + int modelType = m->getModelType(g_ent); + apf::setScalar(modelTagge,ent,0,modelTag); + } + m->end(it); + + apf::Field* modelTyper = apf::createLagrangeField(m,"modelType",apf::SCALAR,1); + it = m->begin(0); + while( (ent = m->iterate(it)) ) + { + apf::ModelEntity* g_ent = m->toModel(ent); + int modelTag = m->getModelTag(g_ent); + int modelType = m->getModelType(g_ent); + apf::setScalar(modelTyper,ent,0,modelType); + } + m->end(it); } @@ -571,6 +683,13 @@ void Enclosure::makeBox3D(gmi_model* model) return; } +void PiercingCylinder::makePiercingCylinder(gmi_model* model) +{ + + sphereRadius = radius; + gmi_add_analytic(model, dim-1, faceID, cylinderFace, faPer, faRan, 0); +} + void Sphere::makeSphere(gmi_model* model) { @@ -579,6 +698,51 @@ void Sphere::makeSphere(gmi_model* model) gmi_add_analytic(model, dim-1, faceID, sphereFace, faPer, faRan, 0); } + +void setParameterizationCylinder(gmi_model* model,apf::Mesh2* m, Enclosure box, PiercingCylinder cylinder,Sphere circle1, Sphere circle2 ) +{ + //Get the classification of each entity in the SimMesh + apf::MeshEntity* ent; + for(int i =0;i<4;i++) + { + apf::MeshIterator* it = m->begin(i); + while( (ent = m->iterate(it))) + { + apf::ModelEntity* g_ent = m->toModel(ent); + int modelTag = m->getModelTag(g_ent); + int modelType = m->getModelType(g_ent); + if(modelTag > 139 && modelTag< 148) + { + m->setModelEntity(ent,(apf::ModelEntity*)gmi_find(model,modelType,cylinder.faceID)); + } + else if(modelTag >= 200){ + //std::cout< 200 && modelTag <= 204) + m->setModelEntity(ent,(apf::ModelEntity*)gmi_find(model,modelType,circle1.faceID)); + else if(modelTag >= 209 && modelTag <= 212) + m->setModelEntity(ent,(apf::ModelEntity*)gmi_find(model,modelType,circle2.faceID)); + } + else{ + //std::cout<<"belongs nowhere? "<setModelEntity(ent,(apf::ModelEntity*)gmi_find(model,modelType,modelTag)); + } + + } + m->end(it); + } + checkEntities(m); + apf::writeVtkFiles("beforeAccept",m); + + m->setModel(model); + m->acceptChanges(); + checkEntities(m); + apf::writeVtkFiles("afterAccept",m); + + + Reparam::reparameterizeEntities(model,m,box,cylinder,circle1,circle2); + +} + //model tags are based off gmsh default outputs... void setParameterization(gmi_model* model,apf::Mesh2* m, Enclosure box, Sphere sphere) { @@ -700,6 +864,93 @@ void Reparam::reparameterizeEntities(gmi_model*model,apf::Mesh2*m,Enclosure box, else reparameterize3D(model,m,box, sphere); } +void Reparam::reparameterizeEntities(gmi_model*model,apf::Mesh2*m,Enclosure box, PiercingCylinder cylinder, Sphere circle1, Sphere circle2){ + //reparameterizeCylinder(model,m,box, cylinder,circle1,circle2); + + //Need to set the parametric coordinates of each of the boundary vertices + std::map edgeParam; + int edgeScales[12] = {0,1,0,1,0,1,0,1,2,2,2,2}; + double edgeLengths[3] = {boxLength,boxWidth,boxHeight}; + for(int i=0;i<12;i++) + { + edgeParam[box.edgeMap[i]] = edgeScales[i]; + } + + std::map faceParam; + int faceScales[6][2] = {{0,2},{1,2},{0,2},{1,2},{0,1},{0,1}}; + for(int i = 0; i<6;i++) + { + faceParam[box.faceMap[i]] = &(faceScales[i]); + } + + apf::MeshIterator* it = m->begin(0); + apf::MeshEntity* ent; + while( (ent = m->iterate(it)) ) + { + apf::ModelEntity* g_ent = m->toModel(ent); + + apf::MeshEntity* ev[2]; + m->getDownward(ent,0,ev); + int modelTag = m->getModelTag(g_ent); + int modelType = m->getModelType(g_ent); + if(modelType<3 && modelType!=0) + { + apf::Vector3 pt; + apf::Vector3 newParam; + m->getPoint(ent,0,pt); + if(modelType==1) + { + if(modelTag == circle1.faceID){ + double argy = (pt[1]-xyz_offset[1]); + double argx = (pt[0]-xyz_offset[0]); + if(argx == 0 && argy ==0) + newParam[0] = 0.0; // not sure if this will happen or if this is right + else + newParam[0] = atan2(argy,argx); + m->setParam(ent,newParam); + } + else if(modelTag == circle2.faceID){ + double argy = (pt[1]-xyz_offset[1]); + double argx = (pt[0]-xyz_offset[0]); + if(argx == 0 && argy ==0) + newParam[0] = 0.0; // not sure if this will happen or if this is right + else + newParam[0] = atan2(argy,argx); + m->setParam(ent,newParam); + } + else{ + int relevantIndex = edgeParam[modelTag]; + newParam[0]=pt[relevantIndex]/edgeLengths[relevantIndex]; + m->setParam(ent,newParam); + } + } + else if (modelType==2 && modelTag!=cylinder.faceID) + { + int* relevantIndex = faceParam[modelTag][0]; //size is 2 + newParam[0] = pt[relevantIndex[0]]/edgeLengths[relevantIndex[0]]; + newParam[1] = pt[relevantIndex[1]]/edgeLengths[relevantIndex[1]]; + m->setParam(ent,newParam); + } + else if (modelType==2 && modelTag == cylinder.faceID) + { + double argy = (pt[1]-xyz_offset[1]); + double argx = (pt[0]-xyz_offset[0]); + if(argx == 0 && argy ==0) + newParam[0] = 0.0; // not sure if this will happen or if this is right + else + newParam[0] = atan2(argy,argx); + newParam[1] = pt[2]/boxHeight; + + m->setParam(ent,newParam); + } + + } //end if + } //end while + m->end(it); + m->acceptChanges(); + +} + void Reparam::reparameterize2D(gmi_model*model,apf::Mesh2*m,Enclosure box, Sphere sphere){ std::map edgeParam; @@ -795,9 +1046,12 @@ void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ //at this point, hmin and hmax haven't been set yet //apf::Field* size_initial = apf::createLagrangeField(m,"size_initial",apf::SCALAR,1); lion_set_verbosity(1); - size_iso = apf::createLagrangeField(m,"proteus_size",apf::SCALAR,1); + apf::Field* size_init = apf::createLagrangeField(m,"proteus_init",apf::SCALAR,1); + sizeFieldList.push(size_init); apf::MeshIterator* it = m->begin(0); apf::MeshEntity* ent; + hmin = 0.0125; + hmax = 0.2; while( (ent = m->iterate(it)) ) { apf::Vector3 pt; @@ -812,36 +1066,26 @@ void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ apf::setScalar(size_iso,ent,0,hmax); } */ - apf::setScalar(size_iso,ent,0,0.25); - } - m->end(it); - - //gradeMesh(1.5); - - apf::Field* modelTagge = apf::createLagrangeField(m,"modelTags",apf::SCALAR,1); - it = m->begin(0); - while( (ent = m->iterate(it)) ) - { - apf::ModelEntity* g_ent = m->toModel(ent); - int modelTag = m->getModelTag(g_ent); - int modelType = m->getModelType(g_ent); - apf::setScalar(modelTagge,ent,0,modelTag); - } - m->end(it); - - apf::Field* modelTyper = apf::createLagrangeField(m,"modelType",apf::SCALAR,1); - it = m->begin(0); - while( (ent = m->iterate(it)) ) - { + //apf::setScalar(size_iso,ent,0,0.2); apf::ModelEntity* g_ent = m->toModel(ent); int modelTag = m->getModelTag(g_ent); int modelType = m->getModelType(g_ent); - apf::setScalar(modelTyper,ent,0,modelType); + if(modelTag == modelPiercingCylinder.faceID && modelType==2 || modelType==1 && (modelTag == modelCircle1.faceID || modelTag == modelCircle2.faceID)) + apf::setScalar(size_init,ent,0,0.0125); + else + apf::setScalar(size_init,ent,0,0.1025); } m->end(it); + apf::writeVtkFiles("pregradeProteus",m); + std::cout<<"grade mesh intiial\n"; + //gradeMesh(1.5); + isotropicIntersect(); + std::cout<<"grade mesh post\n"; + + checkEntities(m); apf::writeVtkFiles("initialProteus",m); - ma::Input* in = ma::configure(m,size_iso); + ma::Input* in = ma::configure(m,size_init); in->maximumIterations = 10; in->shouldSnap = true; in->shouldTransferParametric = true; @@ -852,9 +1096,10 @@ void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ apf::writeVtkFiles("middleProteus",m); //apf::writeVtkFiles("initialAdapt",m); - freeField(size_iso); + freeField(size_init); - size_iso = apf::createLagrangeField(m,"proteus_size",apf::SCALAR,1); + size_init = apf::createLagrangeField(m,"proteus_initial",apf::SCALAR,1); + sizeFieldList.push(size_init); it = m->begin(0); while( (ent = m->iterate(it)) ) { @@ -866,24 +1111,32 @@ void MeshAdaptPUMIDrvr::initialAdapt_analytic(){ else apf::setScalar(size_iso,ent,0,hmax); */ - apf::setScalar(size_iso,ent,0,0.25); + //apf::setScalar(size_iso,ent,0,0.05); + apf::ModelEntity* g_ent = m->toModel(ent); + int modelTag = m->getModelTag(g_ent); + int modelType = m->getModelType(g_ent); + if(modelTag == modelPiercingCylinder.faceID && modelType==2 || modelType==1 && (modelTag == modelCircle1.faceID || modelTag == modelCircle2.faceID)) + apf::setScalar(size_init,ent,0,0.0125); + else + apf::setScalar(size_init,ent,0,0.1025); + } m->end(it); //gradeMesh(1.5); - in = ma::configure(m,size_iso); + in = ma::configure(m,size_init); in->maximumIterations = 10; in->shouldSnap = true; in->shouldTransferParametric = true; in->shouldFixShape = true; - in->debugFolder="./debug_fine"; - ma::adaptVerbose(in,false); + in->debugFolder="./debug_fine2"; + ma::adaptVerbose(in,true); m->verify(); apf::writeVtkFiles("finalProteus",m); //apf::writeVtkFiles("initialAdapt2",m); - freeField(size_iso); + freeField(size_init); } @@ -963,6 +1216,95 @@ void MeshAdaptPUMIDrvr::createAnalyticGeometry(int dim, double* boxDim,double*sp } +void MeshAdaptPUMIDrvr::createAnalyticGeometryCylinder(int dim, double* boxDim,double*sphereCenter, double radius) +{ + boxLength = boxDim[0]; + boxWidth = boxDim[1]; + boxHeight = boxDim[2]; + updateSphereCoordinates(sphereCenter); + + //create analytic model + gmi_model* model = gmi_make_analytic(); + + //add sphere + + PiercingCylinder cylinder = PiercingCylinder(); + cylinder.radius = radius; + + cylinder.makePiercingCylinder(model); + //int edgeID[2] = {200,201}; + int edgeID[2] = {9,10}; + Sphere circle1 = Sphere(2); + circle1.faceID = edgeID[0]; + circle1.radius = radius; + //can't use makeSphere because of paramterization function is not general + gmi_add_analytic(model, 1, circle1.faceID, circleFace0, circle1.faPer, circle1.faRan, 0); + + Sphere circle2 = Sphere(2); + circle2.faceID = edgeID[1]; + circle2.radius = radius; + //circle2.makeSphere(model); + gmi_add_analytic(model, 1, circle2.faceID, circleFace1, circle2.faPer, circle2.faRan, 0); + + //add the box + geomDim = dim; + Enclosure box; + + //add special flag to makeBox to include circles? + //need to classify vertices on these circles independently from the cylindrical face + + box.makeBox3D(model); + //create 2 circles, parameterize it as part of the two faces 1 and 6, faces 4+5 + + //face ids come from faceMap + agm_bdry b; + agm_use regionUse; + + b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(gmi_find(model,2,1))); + regionUse = add_adj(model, b, circle1.faceID); + gmi_add_analytic_reparam(model, regionUse, reparam_Circle0, 0); + + b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(gmi_find(model,2,6))); + regionUse = add_adj(model, b, circle2.faceID); + gmi_add_analytic_reparam(model, regionUse, reparam_Circle1, 0); + + //parameterize the cylinder to the region + b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(gmi_find(model,dim,box.regionID))); + regionUse = add_adj(model, b, cylinder.faceID); + gmi_add_analytic_reparam(model, regionUse, regionFunction, 0); + + b = agm_add_bdry(gmi_analytic_topo(model), agm_from_gmi(gmi_find(model,2,cylinder.faceID))); + regionUse = add_adj(model, b, circle1.faceID); + gmi_add_analytic_reparam(model, regionUse, reparam_CircleCylinder0, 0); + regionUse = add_adj(model, b, circle2.faceID); + gmi_add_analytic_reparam(model, regionUse, reparam_CircleCylinder1, 0); + + setParameterizationCylinder(model,m,box,cylinder,circle1,circle2); + + modelBox = box; + modelCircle1 = circle1; + modelCircle2 = circle2; + modelPiercingCylinder = cylinder; + isAnalytic = 1; + m->verify(); + + apf::Migration* plan = 0; + gmi_register_mesh(); + Splitter::partitionFactor = PCU_Comm_Peers(); + Splitter::switchToOriginals(); + bool isOriginal = ((PCU_Comm_Self() % Splitter::partitionFactor) == 0); + if (isOriginal) { + plan = Splitter::getPlan(m); + } + Splitter::switchToAll(); + m = apf::repeatMdsMesh(m, model, plan, Splitter::partitionFactor); + Parma_PrintPtnStats(m, ""); + + return; + +} + + void Splitter::freeMesh(apf::Mesh* m) { m->destroyNative(); diff --git a/proteus/MeshAdaptPUMI/createAnalyticGeometry.h b/proteus/MeshAdaptPUMI/createAnalyticGeometry.h index 7b0241c8c9..d760cba9a2 100644 --- a/proteus/MeshAdaptPUMI/createAnalyticGeometry.h +++ b/proteus/MeshAdaptPUMI/createAnalyticGeometry.h @@ -71,8 +71,26 @@ class Sphere{ void makeSphere(gmi_model* model); }; +class PiercingCylinder{ + public: + int faceID; //values are dictated by spatial tools + double radius; + double offset[3]; + int dim; + int faPer[2] = {1, 0}; + double faRan[2][2] = {{0,6.28318530718},{0.0,apf::pi}}; + + PiercingCylinder(){ + dim = 3; + faceID=11; + } + void makePiercingCylinder(gmi_model* model); +}; + + namespace Reparam{ void reparameterizeEntities(gmi_model*model,apf::Mesh2*m,Enclosure box, Sphere sphere); + void reparameterizeEntities(gmi_model*model,apf::Mesh2*m,Enclosure box, PiercingCylinder cylinder, Sphere circle1, Sphere circle2); void reparameterize2D(gmi_model*model,apf::Mesh2*m,Enclosure box, Sphere sphere); void reparameterize3D(gmi_model*model,apf::Mesh2*m,Enclosure box, Sphere sphere); }