proteus  1.8.1
C/C++/Fortran libraries
cMeshAdaptPUMI.cpp
Go to the documentation of this file.
1 #include <gmi.h>
2 #include <gmi_mesh.h>
3 #include <gmi_null.h>
4 #include <ma.h>
5 #include <maShape.h>
6 #include <apfMDS.h>
7 #include <PCU.h>
8 #include <apf.h>
9 
10 #include <iostream>
11 #include <fstream>
12 #include <stdio.h>
13 #include <string.h>
14 
15 #include "MeshAdaptPUMI.h"
16 #include <sam.h>
17 #include <samSz.h>
18 
19 #include <apfShape.h>
20 
21 
22 extern double dt_err;
23 #ifdef PROTEUS_USE_SIMMETRIX
24 //PROTEUS_USE_SIMMETRIX is a compiler macro that indicates whether Simmetrix libraries are used
25 //This is defined in proteus/config/default.py and is contingent on the existence of a SIM_INCLUDE_DIR path
26  #include <gmi_sim.h>
27  #include <SimUtil.h>
28  #include <SimModel.h>
29  #include <MeshSim.h>
30  #include <SimMeshTools.h>
31  #define FACE 2
32  pAManager SModel_attManager(pModel model);
33 #endif
34 
40 //MeshAdaptPUMIDrvr::MeshAdaptPUMIDrvr(double Hmax, double Hmin, double HPhi,int AdaptMesh, int NumIter, int NumAdaptSteps,const char* sfConfig, const char* maType,const char* logType, double targetError, double targetElementCount,int reconstructedFlag,double maxAspectRatio, double gradingFact)
41 MeshAdaptPUMIDrvr::MeshAdaptPUMIDrvr()
47 {
48  m = 0;
49  PCU_Comm_Init();
50  PCU_Protect();
51 
52 #ifdef PROTEUS_USE_SIMMETRIX
53  Sim_readLicenseFile(0);
54  SimModel_start();
55  gmi_register_sim();
56 #endif
57  nAdapt=0;
58  nTriggers=0;
59  nEstimate=0;
60  global[0] = global[1] = global[2] = global[3] = 0;
61  local[0] = local[1] = local[2] = local[3] = 0;
62  size_iso = 0;
63  size_scale = 0;
64  size_frame = 0;
65  err_reg = 0;
66  vmsErrH1 = 0;
67  errRho_reg = 0;
68  errRel_reg = 0;
69  error_reference = 0;
70  gmi_register_mesh();
71  gmi_register_null();
72  approximation_order = 2;
73  integration_order = 3;
74  total_error = 0.0;
75  errRho_max = 0.0;
76  rel_err_total = 0.0;
77  exteriorGlobaltoLocalElementBoundariesArray = NULL;
78  modelFileName = NULL;
79  adapt_type_config = "test"; //refers to isotropic or anisotropic, feature effectively turned off for now
80  has_gBC = false;
81  target_element_count = 0;//targetElementCount;
82  domainVolume = 0.0;
83  THRESHOLD = 0.0;
84  initialReconstructed = 0;
85  maxAspect = 2.0;//maxAspectRatio;
86  hasIBM=hasInterface=hasVMS=hasERM=hasAniso=hasAnalyticSphere=useProteus=useProteusAniso=0;
87 }
88 
89 MeshAdaptPUMIDrvr::~MeshAdaptPUMIDrvr()
93 {
94 /*
95  freeField(err_reg);
96  freeField(vmsErrH1);
97  freeField(errRho_reg);
98  freeField(errRel_reg);
99  freeField(error_reference);
100  freeField(size_iso);
101  freeField(size_scale);
102  freeField(size_frame);
103 */
104 /*
105  if(isReconstructed){
106  free(modelVertexMaterial);
107  free(modelBoundaryMaterial);
108  free(modelRegionMaterial);
109  //m->destroyNative();
110  //gmi_destroy(m->getModel());
111  //apf::destroyMesh(m);
112  }
113 */
114  PCU_Comm_Free();
115 #ifdef PROTEUS_USE_SIMMETRIX
116  SimModel_stop();
117  Sim_unregisterAllKeys();
118 #endif
119 }
120 
121 
122 static bool ends_with(std::string const& str, std::string const& ext)
123 {
124  return str.size() >= ext.size() &&
125  str.compare(str.size() - ext.size(), ext.size(), ext) == 0;
126 }
136 int MeshAdaptPUMIDrvr::loadModelAndMesh(const char* modelFile, const char* meshFile)
137 {
138  comm_size = PCU_Comm_Peers();
139  comm_rank = PCU_Comm_Self();
140  if (ends_with(meshFile, ".msh")){
141  m = apf::loadMdsFromGmsh(gmi_load(modelFile), meshFile);
142  std::cout<<"Boundary Condition functionality has not been built in for gmsh yet.\n";
143  }
144  else if (ends_with(modelFile,".smd")){
145  m = apf::loadMdsMesh(modelFile, meshFile);
146  modelFileName=(char *) malloc(sizeof(char) * strlen(modelFile));
147  strcpy(modelFileName,modelFile);
148  getSimmetrixBC();
149  }
150  else{
151  m = apf::loadMdsMesh(modelFile, meshFile);
152  }
153 
154  m->verify();
155  return 0;
156 }
157 
158 int MeshAdaptPUMIDrvr::loadMeshForAnalytic(const char* meshFile,double* boxDim,double* sphereCenter, double radius)
159 {
160  //assume analytic
161  comm_size = PCU_Comm_Peers();
162  comm_rank = PCU_Comm_Self();
163  m = apf::loadMdsMesh(".null", meshFile);
164  m->verify();
165 
166  //create analytic geometry
167  gmi_model* testModel = createSphereInBox(boxDim,sphereCenter,radius);
168  m->verify();
169 
170 
171 /*
172  apf::writeVtkFiles("afterAnalytic",m);
173  std::cout<<"test Model "<<testModel<<" mesh model "<<m->getModel()<<std::endl;
174  std::abort();
175 */
176  return 0;
177 }
178 
179 
180 
181 int MeshAdaptPUMIDrvr::getSimmetrixBC()
189 {
190 
191 #ifdef PROTEUS_USE_SIMMETRIX
192  pGModel model = 0;
193  model=GM_load(modelFileName,NULL,NULL);
194 
195  pAManager attmngr = SModel_attManager(model);
196  pACase acase = AMAN_findCaseByType(attmngr, "problem definition");
197  if (acase){
198  if(comm_rank==0)std::cout<<"Found case, setting the model"<<std::endl;
199  AttCase_setModel(acase,model);
200  has_gBC=true;
201  AttCase_associate(acase,NULL);
202 
203  pGFace gFace;
204  GFIter gfIter = GM_faceIter(model);
205  pAttribute Att[GM_numFaces(model)];
206  int attMap[GM_numFaces(model)];
207  int nF=0;
208 
209  char strAtt[2][25] = {"traction vector","comp3"};
210  int modelEntTag;
211 
212  while(gFace = GFIter_next(gfIter))
213  {
214  if(GEN_attrib((pGEntity)gFace,strAtt[0]))
215  {
216  modelEntTag=GEN_tag((pGEntity)gFace);
217  Att[nF]=GEN_attrib((pGEntity)gFace,strAtt[0]);
218  attMap[nF] = modelEntTag;
219  nF++;
220  }
221  }
222  GFIter_delete(gfIter);
223 
224  apf::MeshIterator* fIter = m->begin(FACE);
225  apf::MeshEntity* fEnt;
226  apf::Vector3 evalPt;
227  int numqpt=0;
228  const int nsd = 3;
229  const int bcFlag[nsd+1] = {0,1,1,1};
230 
231  //assign a label to the BC type tag
232  char label[9],labelflux[4][9],type_flag;
233 
234  sprintf(label,"BCtype");
235  BCtag = m->createIntTag(label,4);
236  for(int idx=0;idx<4;idx++)
237  {
238  if(idx == 0) sprintf(&type_flag,"p");
239  else if(idx == 1) sprintf(&type_flag,"u");
240  else if(idx == 2) sprintf(&type_flag,"v");
241  else if(idx == 3) sprintf(&type_flag,"w");
242  if(idx>0) sprintf(labelflux[idx],"%c_flux",type_flag);
243  }
244 
245  while(fEnt = m->iterate(fIter))
246  {
247  apf::ModelEntity* me=m->toModel(fEnt);
248  modelEntTag = m->getModelTag(me);
249  apf::ModelEntity* boundary_face = m->findModelEntity(FACE,modelEntTag);
250  if(numqpt==0)
251  {
252  apf::MeshElement* testElem = apf::createMeshElement(m,fEnt);
253  numqpt = apf::countIntPoints(testElem,integration_order);
254  for(int idx=1;idx<nsd+1;idx++)
255  fluxtag[idx]= m->createDoubleTag(labelflux[idx],numqpt);
256  apf::destroyMeshElement(testElem);
257  }
258  if(me==boundary_face)
259  {
260  for(int i=0;i<nF;i++)
261  {
262  if(attMap[i]==modelEntTag)
263  {
264  apf::MeshElement* testElem = apf::createMeshElement(m,fEnt);
265  double data[nsd+1][numqpt];
266  for(int k=0; k<numqpt;k++)
267  {
268  apf::getIntPoint(testElem,integration_order,k,evalPt);
269  apf::Vector3 evalPtGlobal;
270  apf::mapLocalToGlobal(testElem,evalPt,evalPtGlobal);
271  double evalPtSim[nsd];
272  evalPtGlobal.toArray(evalPtSim);
273  for(int j=0;j<nsd;j++)
274  data[j+1][k]=AttributeTensor1_evalDS((pAttributeTensor1)Att[i], j,evalPtSim);
275  }
276  m->setIntTag(fEnt,BCtag,&(bcFlag[0]));
277  for(int idx=1;idx<nsd+1;idx++)
278  {
279  m->setDoubleTag(fEnt,fluxtag[idx],data[idx]); //set the quadrature points
280  }
281  apf::destroyMeshElement(testElem);
282  break;
283  } //end if on model
284  else
285  {
286  int dummy[4] = {0,0,0,0};
287  m->setIntTag(fEnt,BCtag,&(dummy[0]));
288  }
289  }//end loop over attributes
290  if(nF==0)
291  {
292  int dummy[4] = {0,0,0,0};
293  m->setIntTag(fEnt,BCtag,&(dummy[0]));
294  }
295  }
296  }//end while
297  m->end(fIter);
298  AMAN_release( attmngr );
299  } else {
300  if(comm_rank==0)
301  std::cout<<"Case not found, no BCs?\n"<<std::endl;
302  //exit(1);
303  }
304 
305  if(comm_rank==0)std::cout<<"Finished reading and storing diffusive flux BCs\n";
306 #endif
307  return 0;
308 }
309 static int countTotal(apf::Mesh* m, int dim)
310 {
311  int total = apf::countOwned(m, dim);
312  PCU_Add_Ints(&total, 1);
313  return total;
314 }
315 
316 
317 int MeshAdaptPUMIDrvr::willErrorAdapt()
318 
326 {
327  int adaptFlag=0;
328  int assertFlag;
329 
330  //get current size field
331  apf::Field* currentField;
332 
333 /*
334  if(!size_iso) //if no previous size field
335  {
336  currentField = samSz::isoSize(m);
337  }
338  else //otherwise use previous size field, remember to reflect this in interfaceAdapt or collapse to a single function
339  {
340  currentField = apf::createFieldOn(m, "currentField", apf::SCALAR);
341  apf::copyData(currentField,size_iso);
342  }
343 */
344  //currentField = apf::createFieldOn(m, "currentField", apf::SCALAR);
345  //apf::copyData(currentField,size_iso);
346  currentField = samSz::isoSize(m);
347 
348  //get error-based size field
349  getERMSizeField(total_error);
350  apf::Field* errorField = sizeFieldList.front();
351  sizeFieldList.pop(); //remove this size field from the queue
352 
353 
354  apf::Field *errorTriggered = apf::createLagrangeField(m, "errorTriggered", apf::SCALAR, 1);
355  //determine if desired mesh is contained in current mesh
356  apf::MeshEntity* ent;
357  apf::MeshIterator* it = m->begin(0);
358  while( (ent = m->iterate(it)) )
359  {
360  double h_current = apf::getScalar(currentField,ent,0);
361  double h_needed = apf::getScalar(errorField,ent,0);
362  if(h_current>h_needed*1.5){
363  adaptFlag+=1;
364  apf::setScalar(errorTriggered,ent,0,h_current/h_needed*1.0);
365  //apf::writeVtkFiles("willErrorAdapt", m);
366  //std::cout<<"What is the ent? "<<localNumber(ent)<<std::endl;
367  //std::exit(1);
368  //break;
369  }
370  else
371  apf::setScalar(errorTriggered,ent,0,-1);
372 
373  }//end while
374  m->end(it);
375 
376  //modify error field to be the ratio
377  it = m->begin(0);
378  while( (ent = m->iterate(it)) )
379  {
380  double h_current = apf::getScalar(currentField,ent,0);
381  double h_needed = apf::getScalar(errorField,ent,0);
382  apf::setScalar(errorField,ent,0,h_current/h_needed*1.0);
383  }//end while
384  m->end(it);
385 
386  assertFlag = adaptFlag;
387  PCU_Add_Ints(&assertFlag,1);
388  //assert(assertFlag ==0 || assertFlag == PCU_Proc_Peers());
389 
390  if(assertFlag>0)
391  {
392  double totalNodes = countTotal(m,0);
393  double triggeredPercentage = assertFlag*100.0/totalNodes;
394  char buffer[50];
395  sprintf(buffer,"Need to error adapt %f%%",triggeredPercentage);
396  logEvent(buffer,3);
397 
398 /*
399  if(nTriggers%10==0)
400  {
401  char namebuffer[50];
402  sprintf(namebuffer,"needErrorAdapt_%i",nTriggers);
403  apf::writeVtkFiles(namebuffer, m);
404  }
405 */
406  nTriggers++;
407  }
408 
409  apf::destroyField(currentField);
410  //apf::destroyField(errorField);
411  apf::destroyField(errorTriggered);
412 
413  return assertFlag;
414 }
415 
416 int MeshAdaptPUMIDrvr::willErrorAdapt_reference()
417 {
418  int adaptFlag=0;
419  int assertFlag;
420 
421  getERMSizeField(total_error);
422  sizeFieldList.pop(); //remove this size field from the queue
423 
424  //if(m->findField("errorTriggered"))
425  // apf::destroyField(m->findField("errorTriggered"));
426  apf::Field* errorTriggered;
427  if(m->findField("errorTriggered"))
428  errorTriggered = m->findField("errorTriggered");
429  else
430  errorTriggered = apf::createField(m,"errorTriggered",apf::SCALAR,apf::getVoronoiShape(m->getDimension(),1));
431 
432  apf::Field* error_current = m->findField("VMSH1");
433  apf::Field* error_reference=NULL;
434 
435 
436  if(m->findField("errorRate"))
437  apf::destroyField(m->findField("errorRate"));
438  apf::Field* errorRateField = apf::createField(m,"errorRate",apf::SCALAR,apf::getVoronoiShape(m->getDimension(),1));
439 
440  //need to set the error trigger field to 0
441  apf::MeshEntity* ent;
442  apf::MeshIterator* it;
443 
444  if(nTriggers == 0)
445  {
446  it = m->begin(m->getDimension());
447  while( (ent = m->iterate(it) ) )
448  {
449  apf::setScalar(errorTriggered,ent,0,-1.0);
450  apf::setScalar(errorRateField,ent,0,0.0);
451  }
452  m->end(it);
453  logEvent("SET ERROR TRIGGERED!",4);
454  }
455 
456 
457  //need to have size field here just to define it
458 
459  //set the error reference field to be the current error field for next time otherwise, retrieve the error_reference field
460  if(!m->findField("error_reference"))
461  {
462  T_reference = T_current;
463  error_reference = apf::createField(m,"error_reference",apf::SCALAR,apf::getVoronoiShape(nsd,1));
464  apf::copyData(error_reference,error_current);
465  logEvent("SUCCESSFULLY COPIED!",4);
466 
467  return 0;
468  }
469  else
470  {
471  error_reference = m->findField("error_reference");
472  }
473 
474  double dt_step = dt_err; //global variable imported from VMS
475 
476  if(m->findField("sizeRatio"))
477  apf::destroyField(m->findField("sizeRatio"));
478  apf::Field* sizeRatioField = apf::createField(m,"sizeRatio",apf::SCALAR,apf::getVoronoiShape(m->getDimension(),1));
479 
480  it = m->begin(m->getDimension());
481 
482  while( (ent = m->iterate(it) ) )
483  {
484  double err_local_current = apf::getScalar(error_current,ent,0);
485  double err_local_ref = apf::getScalar(error_reference,ent,0);
486  double errorRate = (err_local_current-err_local_ref)/(T_current-T_reference);
487  apf::setScalar(errorRateField,ent,0,errorRate);
488  double err_predict = errorRate*delta_T_next + err_local_current;
489 
490  //double sizeRatio_local = apf::getScalar(m->findField("sizeRatio"),ent,0);
491  //if(errorRate > 0 && (err_predict > target_error) && sizeRatio_local>2.0 && nTriggers>=5)
492  //double sizeRatio_local_predict = pow(err_predict/target_error,2.0/(2*1+2));
493 
494  double h_old;
495  double h_new;
496  apf::MeshElement* element = apf::createMeshElement(m, ent);
497 
498  if (m->getDimension() == 2)
499  h_old = apf::computeLargestHeightInTri(m,ent);
500  else
501  h_old = apf::computeShortestHeightInTet(m,ent);
502  h_new = h_old * pow((target_error / err_predict),2.0/(2.0*(1.0)+nsd));
503  //clamp h_new and then compare against h_old
504  if(h_new > hmax)
505  h_new = hmax;
506  if(h_new < hmin)
507  h_new = hmin;
508  double sizeRatio_local_predict = h_old/h_new;
509  apf::setScalar(sizeRatioField,ent,0,sizeRatio_local_predict);
510  apf::destroyMeshElement(element);
511 
512  //apf::setScalar(sizeRatioField,ent,0,size_ratio_local_predict);
513  if(errorRate > 0 && sizeRatio_local_predict>2.0 && nTriggers>= (0.1*numAdaptSteps))
514  {
515  if(apf::getScalar(errorTriggered,ent,0)==1)
516  {
517  adaptFlag = 1;
518  apf::setScalar(errorTriggered,ent,0,2.0);
519  }
520  else
521  apf::setScalar(errorTriggered,ent,0,1.0);
522  }
523  else
524  apf::setScalar(errorTriggered,ent,0,-1.0);
525  }
526  m->end(it);
527 
528  //upper bound on adapt steps
529  if(nTriggers>=2*numAdaptSteps)
530  {
531  logEvent("Adapting because upper limit was reached.",4);
532  adaptFlag = 1;
533  }
534 
535  assertFlag = adaptFlag;
536  PCU_Add_Ints(&assertFlag,1);
537  nTriggers++;
538 
539  //if adapt, modify the error field to be predictive
540  //otherwise, store the current error field to be the next reference field
541  if(assertFlag > 0)
542  {
543  logEvent("Need to error based adapt!!!",4);
544 /*
545  it = m->begin(m->getDimension());
546  while( (ent = m->iterate(it) ) )
547  {
548  double err_local_current = apf::getScalar(error_current,ent,0);
549  double err_local_ref = apf::getScalar(error_reference,ent,0);
550  //double errorRate = (err_local_current-err_local_ref)/(T_current-T_reference);
551  double errorRate = apf::getScalar(errorRateField,ent,0);
552 
553  //double err_predict = errorRate*dt_step*numAdaptSteps + err_local_current;
554  double err_predict = errorRate*dt_step*numAdaptSteps/10.0 + err_local_current;
555  if(err_predict > err_local_current)
556  {
557  apf::setScalar(error_current,ent,0,err_predict);
558  }
559  }
560  m->end(it);
561 */
562  }
563  else
564  {
565  it = m->begin(m->getDimension());
566  while( (ent = m->iterate(it) ) )
567  {
568  double err_local_current = apf::getScalar(error_current,ent,0);
569  //set the reference field for the next step
570  apf::setScalar(error_reference,ent,0,err_local_current);
571  }
572 
573  }
574 
575  T_reference = T_current;
576 
577  return assertFlag;
578 }
579 
580 
581 
582 int MeshAdaptPUMIDrvr::willAdapt()
583 //Master function that calls other adapt-trigger functions
584 {
585  int adaptFlag = 0;
586  //TODO: need to separate IBM from interface to allow for two-phase with IBM
587  if(hasInterface or hasIBM)
588  {
589  adaptFlag += willInterfaceAdapt();
590  logEvent("will adapt because of interface",4);
591  }
592  if(hasVMS)
593  {
594  adaptFlag += willErrorAdapt_reference();
595  logEvent("will adapt because of error",4);
596  }
597  //TODO: need to do AllReduce here instead of in each function
598  if(adaptFlag > 0)
599  adaptFlag = 1;
600 
601  if(adaptFlag == 0)
602  {
603  //allocated in transfer fields... this is not a good way of doing things, but don't know how to pass a numpy array without having to allocate memory just yet
604  free(rho);
605  free(nu);
606  }
607 
608  return adaptFlag;
609 }
610 
611 
612 int MeshAdaptPUMIDrvr::willInterfaceAdapt()
613 //Does banded adapt need to happen for an isotropic mesh?
614 //I need to loop over all mesh edges and determine if the edge intersects the blending region.
615 //If so, need to check the size values on the edge-adjacent vertices.
616 //If either size value is greater than h_interface*1.5, then we know we need to adapt
617 {
618  int adaptFlag=0;
619  int assertFlag;
620 
621  //get current size field
622  apf::Field* currentField;
623 /*
624  if(!size_iso) //if no previous size field
625  {
626  currentField = samSz::isoSize(m);
627  }
628  else //otherwise use previous size field, remember to reflect this in interfaceAdapt or collapse to a single function
629  {
630  currentField = apf::createFieldOn(m, "currentField", apf::SCALAR);
631  apf::copyData(currentField,size_iso);
632  }
633 */
634 
635  currentField = samSz::isoSize(m);
636  double edgeRatio = 1.5; //need to be taken from MeshAdapt library
637 
638  //get banded size field
639  double L_band = (N_interface_band)*hPhi;
640  calculateSizeField(L_band);
641  apf::Field* interfaceField = sizeFieldList.front();
642  sizeFieldList.pop(); //destroy this size field
643 
644 
645  //determine if desired mesh is contained in current mesh
646  apf::MeshEntity* ent;
647  apf::MeshIterator* it = m->begin(0);
648  while( (ent = m->iterate(it)) )
649  {
650  double h_current = apf::getScalar(currentField,ent,0);
651  double h_needed = apf::getScalar(interfaceField,ent,0);
652  if(h_current/h_needed > edgeRatio){
653  adaptFlag=1;
654  break;
655  }
656  }//end while
657 
658  assertFlag = adaptFlag;
659  PCU_Add_Ints(&assertFlag,1);
660  //assert(assertFlag ==0 || assertFlag == PCU_Proc_Peers());
661 
662  apf::destroyField(currentField);
663  apf::destroyField(interfaceField);
664 
665  return assertFlag;
666 }
667 
668 
669 
670 int MeshAdaptPUMIDrvr::adaptPUMIMesh(const char* inputString)
682 {
683  if (hasAniso)
684  calculateAnisoSizeField();
685  if(hasERM)
686  {
687  assert(err_reg);
688  removeBCData();
689  double t1 = PCU_Time();
690  getERMSizeField(total_error);
691  double t2 = PCU_Time();
692  if(comm_rank==0 && logging_config == "on"){
693  std::ofstream myfile;
694  myfile.open("error_estimator_timing.txt", std::ios::app );
695  myfile << t2-t1<<std::endl;
696  myfile.close();
697  }
698  }
699  if((hasVMS) && std::string(inputString)==""){
700  assert(vmsErrH1);
701  getERMSizeField(total_error);
702  }
703  if(hasAnalyticSphere)
704  {
705  setSphereSizeField();
706  }
707  if(hasInterface || hasIBM)
708  {
709  double L_band = (N_interface_band+1)*hPhi;
710  calculateSizeField(L_band);
711  if(nAdapt>1 && std::string(inputString)=="interface")
712  predictiveInterfacePropagation();
713  else if(nAdapt > 2 && std::string(inputString)=="")
714  predictiveInterfacePropagation();
715  }
716  if (useProteus)
717  size_iso = m->findField("proteus_size");
718  if (useProteusAniso){
719  size_frame = m->findField("proteus_sizeFrame");
720  size_scale = m->findField("proteus_sizeScale");
721  adapt_type_config = "anisotropic";
722  }
723 /*
724  if (hasTest)
725  testIsotropicSizeField();
726 */
727 /*
728  if(size_field_config == "uniform"){
729  //special situation where I only care about err_reg
730  freeField(errRho_reg);
731  freeField(errRel_reg);
732  }
733 */
734  isotropicIntersect();
735 
736  if(logging_config=="on"){
737  char namebuffer[50];
738  sprintf(namebuffer,"pumi_preadapt_%i",nAdapt);
739  apf::writeVtkFiles(namebuffer, m);
740  //sprintf(namebuffer,"beforeAnisotropicAdapt%i_.smb",nAdapt);
741  //m->writeNative(namebuffer);
742  }
743 
744  if(hasERM){
745  //MeshAdapt error will be thrown if region fields are not freed
746  freeField(err_reg);
747  freeField(errRho_reg);
748  freeField(errRel_reg);
749  }
750  if(hasVMS){
751  freeField(vmsErrH1);
752  if(PCU_Comm_Self()==0) std::cout<<"cleared VMS field\n";
753  }
754 
755  // These are relics from an attempt to pass BCs from proteus into the error estimator.
756  // They maybe useful in the future.
757  //m->destroyTag(fluxtag[1]); m->destroyTag(fluxtag[2]); m->destroyTag(fluxtag[3]);
758  delete [] exteriorGlobaltoLocalElementBoundariesArray;
759  exteriorGlobaltoLocalElementBoundariesArray = NULL;
760 
761  for (int d = 0; d <= m->getDimension(); ++d)
762  freeNumbering(local[d]);
763 
764  apf::Field* adaptSize;
765  apf::Field* adaptFrame;
766 
768  ma::Input* in;
769  if(size_field_config == "uniform"){
770  in = ma::makeAdvanced(ma::configureUniformRefine(m));
771  in->shouldFixShape=false;
772  }
773  else{
774  assert(size_iso || (size_scale && size_frame));
775  if(adapt_type_config=="anisotropic" || size_field_config== "interface"){
776  //in = ma::configure(m, size_scale, size_frame);
777  adaptSize = apf::createFieldOn(m, "adapt_size", apf::VECTOR);
778  adaptFrame = apf::createFieldOn(m, "adapt_frame", apf::MATRIX);
779  apf::copyData(adaptSize, size_scale);
780  apf::copyData(adaptFrame, size_frame);
781  in = ma::makeAdvanced(ma::configure(m, adaptSize, adaptFrame));
782  }
783  else{
784  adaptSize = apf::createFieldOn(m, "adapt_size", apf::SCALAR);
785  apf::copyData(adaptSize, size_iso);
786  in = ma::makeAdvanced(ma::configure(m, adaptSize));
787  }
788  }
789 
790  ma::validateInput(in);
791  in->shouldRunPreZoltan = true;
792  in->shouldRunMidZoltan = true;
793  in->shouldRunPostZoltan = true;
794  in->maximumImbalance = 1.05;
795  in->maximumIterations = numIter;
796  if(size_field_config == "meshQuality")
797  {
798  in->shouldSnap = true;
799  in->shouldTransferParametric=true;
800  }
801  else
802  in->shouldSnap = false;
803  //in->goodQuality = 0.16;//0.027;
804  //double mass_before = getTotalMass();
805 
806  double t1 = PCU_Time();
807  //ma::adapt(in);
808  ma::adaptVerbose(in);
809  double t2 = PCU_Time();
810 
811  m->verify();
812  //double mass_after = getTotalMass();
813  //PCU_Add_Doubles(&mass_before,1);
814  //PCU_Add_Doubles(&mass_after,1);
815  if(comm_rank==0 && logging_config=="on"){
816 /*
817  std::ios::fmtflags saved(std::cout.flags());
818  std::cout<<std::setprecision(15)<<"Mass Before "<<mass_before<<" After "<<mass_after<<" diff "<<mass_after-mass_before<<std::endl;
819  std::cout.flags(saved);
820  std::ofstream myfile;
821  myfile.open("adapt_timing.txt", std::ios::app);
822  myfile << t2-t1<<std::endl;
823  myfile.close();
824  std::ofstream mymass;
825  mymass.open("mass_check.txt", std::ios::app);
826  mymass <<std::setprecision(15)<<mass_before<<","<<mass_after<<","<<mass_after-mass_before<<std::endl;
827  mymass.close();
828 */
829  }
830 
831  if(hasERM){
832  if (has_gBC)
833  getSimmetrixBC();
834  }
835  if(logging_config=="on"){
836  char namebuffer[50];
837  sprintf(namebuffer,"pumi_postadapt_%i",nAdapt);
838  apf::writeVtkFiles(namebuffer, m);
839  //sprintf(namebuffer,"afterAnisotropicAdapt%i_.smb",nAdapt);
840  //m->writeNative(namebuffer);
841  }
842  //isReconstructed = 0; //this is needed to maintain consistency with the post-adapt conversion back to Proteus
843  apf::destroyField(adaptSize);
844  if(adapt_type_config=="anisotropic")
845  apf::destroyField(adaptFrame);
846  nAdapt++; //counter for number of adapt steps
847 
848  if(logging_config=="debugRestart")
849  m->writeNative("DEBUG_restart.smb");
850 
851  nTriggers=0;
852  return 0;
853 }
854 
855 double MeshAdaptPUMIDrvr::getMinimumQuality()
862 {
863  ma::SizeField* isf = new ma::IdentitySizeField(m);
864  apf::MeshIterator* it = m->begin(m->getDimension());
865  apf::MeshEntity* e;
866  double minq = 1;
867  while ((e = m->iterate(it)))
868  minq = std::min(minq, ma::measureElementQuality(m, isf, e));
869  m->end(it);
870  delete isf;
871  return PCU_Min_Double(minq);
872 }
873 
874 double MeshAdaptPUMIDrvr::getTotalMass()
880 {
881  apf::Field* voff = m->findField("vof");
882  assert(voff);
883  apf::MeshEntity* e;
884  apf::MeshIterator* it = m->begin(m->getDimension());
885  double mass = 0.0;
886  while ((e = m->iterate(it))) {
887  apf::MeshElement* elem = apf::createMeshElement(m,e);
888  apf::Element* voff_elem = apf::createElement(voff, elem);
889  int int_order = 4;
890  for(int l = 0; l < apf::countIntPoints(elem, int_order); ++l) {
891  apf::Vector3 qpt;
892  apf::getIntPoint(elem,int_order,l,qpt);
893  double vof_val = apf::getScalar(voff_elem,qpt);
894  double rho_val = getMPvalue(vof_val,rho[0],rho[1]);
895  double weight = apf::getIntWeight(elem,int_order,l);
896  apf::Matrix3x3 J;
897  apf::getJacobian(elem,qpt,J); //evaluate the Jacobian at the quadrature point
898  double Jdet = apf::getJacobianDeterminant(J,m->getDimension());
899  mass += rho_val*weight*Jdet;
900  }
901  apf::destroyElement(voff_elem);
902  apf::destroyMeshElement(elem);
903  }
904  m->end(it);
905  PCU_Add_Doubles(&mass,1);
906  return mass;
907 }
910 //Save mesh with solution
911 
912 void MeshAdaptPUMIDrvr::writeMesh(const char* meshFile)
913 {
914  m->writeNative(meshFile);
915  //apf::writeVtkFiles(meshFile,m);
916 }
917 
918 //Clean mesh of all fields and tags
919 
920 void MeshAdaptPUMIDrvr::cleanMesh()
921 {
922  //destroy all fields...
923 
924  for(int i =0;i<m->countFields();i++)
925  {
926  apf::Field* sample = m->getField(i);
927  freeField(sample);
928  }
929  //std::cout<<"find field " <<m->getField(m->countFields())<<" how many fields? "<<m->countFields()<<std::endl;
930  //std::cout<<"is velocity_old here? "<<m->findField("velocity_old")<<std::endl;
931  apf::Field* sample = m->findField("velocity_old");
932  freeField(sample);
933 
934  sample = m->findField("vof_old");
935  freeField(sample);
936  sample = m->findField("ls_old");
937  freeField(sample);
938  sample = m->findField("phi");
939  freeField(sample);
940  sample = m->findField("phi_old");
941  freeField(sample);
942  sample = m->findField("phi_old_old");
943  freeField(sample);
944  sample = m->findField("phid_old");
945  freeField(sample);
946  sample = m->findField("phiCorr");
947  freeField(sample);
948  sample = m->findField("phiCorr_old");
949  freeField(sample);
950  sample = m->findField("phiCorr_old_old");
951  freeField(sample);
952  sample = m->findField("p_old");
953  freeField(sample);
954  sample = m->findField("p");
955  freeField(sample);
956  sample = m->findField("p_old_old");
957  freeField(sample);
958 
959  sample = m->findField("VMSH1");
960  freeField(sample);
961  sample = m->findField("VMSL2");
962  freeField(sample);
963 
964  //destroy all tags
965  apf::DynamicArray<apf::MeshTag*> listTags;
966  m->getTags(listTags);
967  int nTags = listTags.getSize();
968  int numDim = m->getDimension();
969  for(int i=0; i < nTags; i++)
970  {
971  std::string ignoreString ("proteus_number");
972  std::string tagName (m->getTagName(listTags[i]));
973  if(tagName.find(ignoreString) != std::string::npos)
974  {
975  //do nothing
976  }
977  else{
978  for(int j=0;j<(numDim+1);j++)
979  {
980  apf::MeshIterator* it = m->begin(j);
981  apf::MeshEntity* ent;
982  while( (ent = m->iterate(it)) )
983  {
984  if(m->hasTag(ent,listTags[i]))
985  m->removeTag(ent,listTags[i]);
986  }
987  m->end(it);
988  }
989  m->destroyTag(listTags[i]);
990  }
991  }
992 }
993 
994 void MeshAdaptPUMIDrvr::set_nAdapt(int numberAdapt)
995 {
996  nAdapt = numberAdapt;
997  return;
998 }
999 
1000 int MeshAdaptPUMIDrvr::setAdaptProperties(std::vector<std::string> sizeInputs,bool in_adapt, double in_hmax,double in_hmin,double in_hphi, int in_numAdaptSteps, double in_targetError, double in_gradingFactor, bool in_logging, int in_numIterations)
1001 {
1002  for (std::vector<std::string>::iterator it = sizeInputs.begin() ; it != sizeInputs.end(); ++it)
1003  {
1004  if(*it == "ibm")
1005  hasIBM=1;
1006  if(*it == "interface")
1007  hasInterface=1;
1008  if(*it == "error_vms")
1009  hasVMS=1;
1010  if(*it == "error_erm")
1011  hasERM=1;
1012  }
1013  hmin=in_hmin;
1014  hmax=in_hmax;
1015  hPhi=in_hphi;
1016  numIter=in_numIterations;
1017  adaptMesh = in_adapt;
1018  numAdaptSteps = in_numAdaptSteps;
1019  if(PCU_Comm_Self()==0)
1020  printf("MeshAdapt: Setting hmax=%lf, hmin=%lf, numIters(meshadapt)=%d\n",
1021  hmax, hmin, numIter);
1022  logging_config = in_logging;
1023  target_error = in_targetError;
1024  gradingFactor = in_gradingFactor;
1025 
1026  return 0;
1027 }
MeshAdaptPUMIDrvr::~MeshAdaptPUMIDrvr
~MeshAdaptPUMIDrvr()
Definition: cMeshAdaptPUMI.cpp:89
MeshAdaptPUMIDrvr::hasERM
bool hasERM
Definition: MeshAdaptPUMI.h:105
MeshAdaptPUMIDrvr::hPhi
double hPhi
Definition: MeshAdaptPUMI.h:91
MeshAdaptPUMIDrvr::hmin
double hmin
Definition: MeshAdaptPUMI.h:91
MeshAdaptPUMIDrvr::total_error
double total_error
Definition: MeshAdaptPUMI.h:141
MeshAdaptPUMI.h
MeshAdaptPUMIDrvr::MeshAdaptPUMIDrvr
MeshAdaptPUMIDrvr()
Definition: cMeshAdaptPUMI.cpp:41
MeshAdaptPUMIDrvr::predictiveInterfacePropagation
void predictiveInterfacePropagation()
Definition: SizeField.cpp:297
MeshAdaptPUMIDrvr::setSphereSizeField
int setSphereSizeField()
Definition: SizeField.cpp:61
MeshAdaptPUMIDrvr::nsd
int nsd
Definition: MeshAdaptPUMI.h:96
MeshAdaptPUMIDrvr::willErrorAdapt_reference
int willErrorAdapt_reference()
Definition: cMeshAdaptPUMI.cpp:416
MeshAdaptPUMIDrvr::getSimmetrixBC
int getSimmetrixBC()
Function used to read in diffusive flux BC from Simmetrix Model.
Definition: cMeshAdaptPUMI.cpp:181
MeshAdaptPUMIDrvr::hmax
double hmax
Definition: MeshAdaptPUMI.h:91
MeshAdaptPUMIDrvr::calculateAnisoSizeField
int calculateAnisoSizeField()
Definition: SizeField.cpp:1055
MeshAdaptPUMIDrvr::initialReconstructed
int initialReconstructed
Definition: MeshAdaptPUMI.h:147
MeshAdaptPUMIDrvr::hasIBM
bool hasIBM
Definition: MeshAdaptPUMI.h:102
MeshAdaptPUMIDrvr::hasAniso
bool hasAniso
Definition: MeshAdaptPUMI.h:106
MeshAdaptPUMIDrvr::getMPvalue
double getMPvalue(double field_val, double val_0, double val_1)
Function primarily used to get the VOF-weighted average of physical properties at a given point.
Definition: ErrorResidualMethod.cpp:35
MeshAdaptPUMIDrvr::getERMSizeField
int getERMSizeField(double err_total)
Definition: SizeField.cpp:1195
MeshAdaptPUMIDrvr::hasAnalyticSphere
bool hasAnalyticSphere
Definition: MeshAdaptPUMI.h:107
MeshAdaptPUMIDrvr::createSphereInBox
gmi_model * createSphereInBox(double *boxDim, double *sphereCenter, double radius)
Definition: createAnalyticGeometry.cpp:602
MeshAdaptPUMIDrvr::getMinimumQuality
double getMinimumQuality()
Function used to get the worst element quality in the mesh.
Definition: cMeshAdaptPUMI.cpp:855
MeshAdaptPUMIDrvr::writeMesh
void writeMesh(const char *meshFile)
Definition: cMeshAdaptPUMI.cpp:912
MeshAdaptPUMIDrvr::hasVMS
bool hasVMS
Definition: MeshAdaptPUMI.h:104
MeshAdaptPUMIDrvr::maxAspect
int maxAspect
Definition: MeshAdaptPUMI.h:97
MeshAdaptPUMIDrvr::getTotalMass
double getTotalMass()
Function to track total mass of the domain.
Definition: cMeshAdaptPUMI.cpp:874
MeshAdaptPUMIDrvr::rel_err_total
double rel_err_total
Definition: MeshAdaptPUMI.h:143
MeshAdaptPUMIDrvr::nEstimate
int nEstimate
Definition: MeshAdaptPUMI.h:95
MeshAdaptPUMIDrvr::N_interface_band
double N_interface_band
Definition: MeshAdaptPUMI.h:100
MeshAdaptPUMIDrvr::exteriorGlobaltoLocalElementBoundariesArray
int * exteriorGlobaltoLocalElementBoundariesArray
Definition: MeshAdaptPUMI.h:134
MeshAdaptPUMIDrvr::set_nAdapt
void set_nAdapt(int numberAdapt)
Definition: cMeshAdaptPUMI.cpp:994
min
#define min(a, b)
Definition: jf.h:71
MeshAdaptPUMIDrvr::loadModelAndMesh
int loadModelAndMesh(const char *modelFile, const char *meshFile)
Load the mesh and model for SCOREC libraries.
Definition: cMeshAdaptPUMI.cpp:136
dt_err
double dt_err
Definition: VMS.cpp:18
MeshAdaptPUMIDrvr::gradingFactor
double gradingFactor
Definition: MeshAdaptPUMI.h:101
MeshAdaptPUMIDrvr::loadMeshForAnalytic
int loadMeshForAnalytic(const char *meshFile, double *boxDim, double *sphereCenter, double radius)
Definition: cMeshAdaptPUMI.cpp:158
MeshAdaptPUMIDrvr::nTriggers
int nTriggers
Definition: MeshAdaptPUMI.h:94
MeshAdaptPUMIDrvr::integration_order
int integration_order
Definition: MeshAdaptPUMI.h:138
MeshAdaptPUMIDrvr::adaptPUMIMesh
int adaptPUMIMesh(const char *input)
Function used to trigger adaptation.
Definition: cMeshAdaptPUMI.cpp:670
MeshAdaptPUMIDrvr::approximation_order
int approximation_order
Definition: MeshAdaptPUMI.h:137
MeshAdaptPUMIDrvr::useProteusAniso
bool useProteusAniso
Definition: MeshAdaptPUMI.h:109
MeshAdaptPUMIDrvr::modelFileName
char * modelFileName
Definition: MeshAdaptPUMI.h:124
MeshAdaptPUMIDrvr::willErrorAdapt
int willErrorAdapt()
Looks at the estimated error and determines if mesh adaptation is necessary.
Definition: cMeshAdaptPUMI.cpp:317
MeshAdaptPUMIDrvr::hasInterface
bool hasInterface
Definition: MeshAdaptPUMI.h:103
MeshAdaptPUMIDrvr::errRho_max
double errRho_max
Definition: MeshAdaptPUMI.h:142
int_order
int int_order
Definition: ErrorResidualMethod.cpp:21
MeshAdaptPUMIDrvr::numIter
int numIter
Definition: MeshAdaptPUMI.h:92
MeshAdaptPUMIDrvr::size_field_config
std::string size_field_config
Definition: MeshAdaptPUMI.h:114
MeshAdaptPUMIDrvr::willInterfaceAdapt
int willInterfaceAdapt()
Definition: cMeshAdaptPUMI.cpp:612
MeshAdaptPUMIDrvr::BCtag
apf::MeshTag * BCtag
Definition: MeshAdaptPUMI.h:130
MeshAdaptPUMIDrvr::useProteus
bool useProteus
Definition: MeshAdaptPUMI.h:108
MeshAdaptPUMIDrvr::numAdaptSteps
int numAdaptSteps
Definition: MeshAdaptPUMI.h:99
MeshAdaptPUMIDrvr::setAdaptProperties
int setAdaptProperties(std::vector< std::string > sizeInputs, bool in_adapt, double in_hmax, double in_hmin, double in_hphi, int in_numAdaptSteps, double in_targetError, double in_gradingFactor, bool in_logging, int in_numIterations)
Definition: cMeshAdaptPUMI.cpp:1000
MeshAdaptPUMIDrvr::adaptMesh
int adaptMesh
Definition: MeshAdaptPUMI.h:98
MeshAdaptPUMIDrvr::calculateSizeField
int calculateSizeField(double L_band)
Definition: SizeField.cpp:84
MeshAdaptPUMIDrvr::nAdapt
int nAdapt
Definition: MeshAdaptPUMI.h:93
MeshAdaptPUMIDrvr::willAdapt
int willAdapt()
Definition: cMeshAdaptPUMI.cpp:582
MeshAdaptPUMIDrvr::fluxtag
apf::MeshTag * fluxtag[4]
Definition: MeshAdaptPUMI.h:132
MeshAdaptPUMIDrvr::cleanMesh
void cleanMesh()
Definition: cMeshAdaptPUMI.cpp:920
MeshAdaptPUMIDrvr::logging_config
std::string logging_config
Definition: MeshAdaptPUMI.h:116
MeshAdaptPUMIDrvr::adapt_type_config
std::string adapt_type_config
Definition: MeshAdaptPUMI.h:115
MeshAdaptPUMIDrvr::removeBCData
void removeBCData()
Function used to remove the BC tags that were created during the computeDiffusiveFlux() function.
Definition: ErrorResidualMethod.cpp:626
Generated on Fri Jul 1 2022 11:15:47 for proteus by doxygen 1.8.20