2D Analytic Geometry using PUMI

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,19 @@ 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'] != []:
+                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)
@@ -662,20 +709,14 @@ 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")
 
-              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:

proteus/MeshAdaptPUMI/ErrorResidualMethod.cpp

@@ -944,7 +944,7 @@ void MeshAdaptPUMIDrvr::get_local_error(double &total_error)
     std::cerr<<"Error density maximum "<<errRho_max<<std::endl;
     std::cerr<<"U_norm_total "<<u_norm_total<<std::endl;
   }
-
+/*
   if(logging_config=="errorOnly"){ //feature to just look at the error fields without adapting the mesh
     if(comm_rank==0)
       std::cout<<"outputting error field\n";
@@ -956,6 +956,7 @@ void MeshAdaptPUMIDrvr::get_local_error(double &total_error)
     nEstimate++;
     removeBCData();
   }
+*/
 
   m->end(iter);
   apf::destroyField(visc);

proteus/MeshAdaptPUMI/MeshAdaptPUMI.h

@@ -4,6 +4,7 @@
 #include <apfNumbering.h>
 #include <queue>
 #include "PyEmbeddedFunctions.h"
+#include "createAnalyticGeometry.h"
 
 /**
    \file MeshAdaptPUMI.h
@@ -73,8 +74,16 @@ class MeshAdaptPUMIDrvr{
   int gradeMesh(double gradationFactor);
 
   //analytic geometry
-  gmi_model* createSphereInBox(double* boxDim, double*sphereCenter,double radius);
+  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();
 
   //Quality Check Functions
   double getMinimumQuality();
@@ -107,13 +116,16 @@ class MeshAdaptPUMIDrvr{
   bool hasAnalyticSphere;
   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);

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;
 }
 

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);

proteus/MeshAdaptPUMI/cMeshAdaptPUMI.cpp

@@ -84,6 +84,8 @@ MeshAdaptPUMIDrvr::MeshAdaptPUMIDrvr()
   initialReconstructed = 0;
   maxAspect = 2.0;//maxAspectRatio;
   hasIBM=hasInterface=hasVMS=hasERM=hasAniso=hasAnalyticSphere=useProteus=useProteusAniso=0;
+  isAnalytic=0;
+  useQuality=0;
 }
 
 MeshAdaptPUMIDrvr::~MeshAdaptPUMIDrvr()
@@ -160,19 +162,18 @@ 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 = createSphereInBox(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();
 
-
-/*
-  apf::writeVtkFiles("afterAnalytic",m);
-  std::cout<<"test Model "<<testModel<<" mesh model "<<m->getModel()<<std::endl;
-  std::abort();
-*/
   return 0;
 }
 
@@ -689,7 +690,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<<std::endl;
@@ -720,20 +721,33 @@ int MeshAdaptPUMIDrvr::adaptPUMIMesh(const char* inputString)
     size_scale = m->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);
+    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();
 
-  if(logging_config=="on"){
+  if(logging_config){
     char namebuffer[50];
     sprintf(namebuffer,"pumi_preadapt_%i",nAdapt);
     apf::writeVtkFiles(namebuffer, m);
@@ -793,26 +807,24 @@ 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=0;
+  //in->shouldTransferParametric=0;
+  in->shouldSnap=isAnalytic;
+  in->shouldTransferParametric=isAnalytic;
   //in->goodQuality = 0.16;//0.027;
   //double mass_before = getTotalMass();
 
   double t1 = PCU_Time();
   //ma::adapt(in);
-  ma::adaptVerbose(in);
+  in->debugFolder="./debug_fine";
+  ma::adaptVerbose(in,true);
   double t2 = PCU_Time();
 
   m->verify();
   //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<<std::setprecision(15)<<"Mass Before "<<mass_before<<" After "<<mass_after<<" diff "<<mass_after-mass_before<<std::endl;
@@ -832,7 +844,7 @@ int MeshAdaptPUMIDrvr::adaptPUMIMesh(const char* inputString)
     if (has_gBC)
       getSimmetrixBC();
   }
-  if(logging_config=="on"){
+  if(logging_config){
     char namebuffer[50];
     sprintf(namebuffer,"pumi_postadapt_%i",nAdapt);
     apf::writeVtkFiles(namebuffer, m);
@@ -845,10 +857,13 @@ int MeshAdaptPUMIDrvr::adaptPUMIMesh(const char* inputString)
     apf::destroyField(adaptFrame);
   nAdapt++; //counter for number of adapt steps
 
+/*
   if(logging_config=="debugRestart")
     m->writeNative("DEBUG_restart.smb");
+*/
 
   nTriggers=0;
+  //std::exit(1);
   return 0;
 }
 
@@ -1009,6 +1024,8 @@ int MeshAdaptPUMIDrvr::setAdaptProperties(std::vector<std::string> sizeInputs,bo
             hasVMS=1;
         if(*it == "error_erm")
             hasERM=1;
+        if(*it == "meshQuality")
+            useQuality=1;
     }
     hmin=in_hmin; 
     hmax=in_hmax; 
@@ -1025,3 +1042,5 @@ int MeshAdaptPUMIDrvr::setAdaptProperties(std::vector<std::string> sizeInputs,bo
 
     return 0;
 }
+
+