Updated FMM

src/domain/domain.cpp

@@ -139,9 +139,11 @@ std::shared_ptr<ablate::domain::SubDomain> ablate::domain::Domain::GetSubDomain(
         return subDomains.front();
     }
 }
-
+#include "levelSet/interfaceReconstruction.hpp"
+#include "finiteVolume/finiteVolumeSolver.hpp"
 void ablate::domain::Domain::InitializeSubDomains(const std::vector<std::shared_ptr<solver::Solver>>& solvers, const std::shared_ptr<ablate::domain::Initializer>& initializations,
                                                   const std::vector<std::shared_ptr<mathFunctions::FieldFunction>>& exactSolutions) {
+
     // determine the number of fields
     for (auto& solver : solvers) {
         solver->Register(GetSubDomain(solver->GetRegion()));
@@ -158,6 +160,22 @@ void ablate::domain::Domain::InitializeSubDomains(const std::vector<std::shared_
     CreateStructures();
     for (auto& subDomain : subDomains) {
         subDomain->CreateSubDomainStructures();
+
+
+
+auto flow = std::dynamic_pointer_cast<ablate::finiteVolume::FiniteVolumeSolver>(solvers[0]);
+
+
+std::shared_ptr<ablate::levelSet::Reconstruction> reconstruction = std::make_shared<ablate::levelSet::Reconstruction>(subDomain, flow->GetRegionWithoutGhost());
+
+Vec auxVector = subDomain->GetAuxVector(); // Getting the Auxiliary Vector (contains data for all aux fields)
+DM aux_dm = subDomain->GetAuxDM(); // Getting the Auxiliary DM which has Auxiliary fields
+const ablate::domain::Field *levelSetField = &(subDomain->GetField("levelSet")); // Getting the level set field for vertices (FEM)
+reconstruction->arbit_interface(aux_dm, *levelSetField, auxVector);
+
+printf("%s::%d\n", __FILE__, __LINE__);
+exit(0);
+
     }
 
     // set all values to nan to allow for a output check

src/finiteVolume/processes/surfaceForce.cpp

@@ -37,6 +37,14 @@ void ablate::finiteVolume::processes::SurfaceForce::Initialize(ablate::finiteVol
   }
 
   SurfaceForce::reconstruction = std::make_shared<ablate::levelSet::Reconstruction>(SurfaceForce::subDomain, flow.GetRegionWithoutGhost());
+
+
+
+//Vec auxVector = subDomain->GetAuxVector(); // Getting the Auxiliary Vector (contains data for all aux fields)
+//DM aux_dm = subDomain->GetAuxDM(); // Getting the Auxiliary DM which has Auxiliary fields
+//const ablate::domain::Field *levelSetField = &(subDomain->GetField("levelSet")); // Getting the level set field for vertices (FEM)
+//process->reconstruction->arbit_interface(aux_dm, *levelSetField, auxVector);
+
 }
 
 
@@ -121,8 +129,13 @@ PetscErrorCode ablate::finiteVolume::processes::SurfaceForce::ComputeSource(cons
     ablate::domain::Range cellRange;
 
 
-const ablate::domain::Field *vofField = &(subDomain->GetField(TwoPhaseEulerAdvection::VOLUME_FRACTION_FIELD));
-process->reconstruction->ToLevelSet(dm, locX, *vofField);
+//const ablate::domain::Field *vofField = &(subDomain->GetField(TwoPhaseEulerAdvection::VOLUME_FRACTION_FIELD));
+//process->reconstruction->ToLevelSet(dm, locX, *vofField);
+xexit("");
+  Vec auxVector = subDomain->GetAuxVector(); // Getting the Auxiliary Vector (contains data for all aux fields)
+  DM aux_dm = subDomain->GetAuxDM(); // Getting the Auxiliary DM which has Auxiliary fields
+  const ablate::domain::Field *levelSetField = &(subDomain->GetField("levelSet")); // Getting the level set field for vertices (FEM)
+  process->reconstruction->arbit_interface(aux_dm, *levelSetField, auxVector);
 
     const ablate::domain::Field *eulerField = &(subDomain->GetField(ablate::finiteVolume::CompressibleFlowFields::EULER_FIELD));
     DM eulerDM = subDomain->GetFieldDM(*eulerField); // Get an euler-specific DM in case it's not in the same solution vector as the VOF field

src/levelSet/interfaceReconstruction.hpp

@@ -25,7 +25,7 @@ namespace ablate::levelSet {
 
     private:
 
-      const PetscInt nLevels = 12;
+      const PetscInt nLevels = 4000;
 
       enum VecLoc { LOCAL , GLOBAL };
 
@@ -110,14 +110,18 @@ namespace ablate::levelSet {
       void Smooth(const PetscInt *cellMask, const PetscInt *vertMask, Vec lsVec[2], Vec fVec[2]);
 
       void FMM(const PetscInt *cellMask, const PetscInt *vertMask, Vec lsVec[2]);
-      void FMM_CellBased(const PetscInt currentLevel, const PetscInt *cellMask, const PetscInt *vertMask, Vec updatedVec[2], Vec lsVec[2]);
+      void FMM_Hybrid(const PetscInt currentLevel, const PetscInt minVerts, const PetscInt *cellMask, const PetscInt *vertMask, Vec updatedVec[2], Vec lsVec[2], Vec lsVecCopy[2], FILE *animFile) ;
+      void FMM_PrimeHybrid(const PetscInt currentLevel, const PetscInt minVerts, const PetscInt *cellMask, const PetscInt *vertMask, Vec updatedVec[2], Vec lsVec[2], Vec lsVecCopy[2], Vec xVec[2], Vec yVec[2], FILE *animFile) ;
+      void FMM_CellBased(const PetscInt currentLevel, const PetscInt *cellMask, const PetscInt *vertMask, Vec updatedVec[2], Vec lsVec[2], Vec lsVecCopy[2], FILE *animFile);
       void FMM_CellBased_V2(const PetscInt currentLevel, const PetscInt *cellMask, const PetscInt *vertMask, Vec updatedVec[2], Vec lsVec[2]);
       void FMM_CellBased_V3(const PetscInt currentLevel, const PetscInt *cellMask, const PetscInt *vertMask, Vec updatedVec[2], Vec lsVec[2]);
-      void FMM_VertexBased_V1(const PetscInt currentLevel, const PetscInt minVerts, const PetscInt *cellMask, const PetscInt *vertMask, Vec updatedVec[2], Vec lsVec[2]);
-      void FMM_VertexBased_V2(const PetscInt currentLevel, const PetscInt *cellMask, const PetscInt *vertMask, Vec updatedVec[2], Vec lsVec[2]);
-      PetscInt FFM_VertexBased_Solve(const PetscInt dim, const PetscInt minVerts, const PetscReal x0[], const PetscInt nVert, PetscInt verts[], PetscScalar *updatedVertex, PetscScalar *lsArray, PetscReal *updatedLS);
-
-
+      void FMM_VertexBased_V1(const PetscInt currentLevel, const PetscInt minVerts, const PetscInt *cellMask, const PetscInt *vertMask, Vec updatedVec[2], Vec lsVec[2], Vec lsVecCopy[2], FILE *animFile);
+      void FMM_VertexBased_V2(const PetscInt currentLevel, const PetscInt *cellMask, const PetscInt *vertMask, Vec updatedVec[2], Vec lsVec[2], Vec lsVecCopy[2]);
+      PetscInt FFM_VertexBased_Solve(const PetscInt dim, const PetscInt minVerts, const PetscReal x0[], const PetscInt nVert, PetscInt verts[], PetscScalar *updatedVertex, PetscScalar *lsArray, PetscReal *updatedLS, const PetscInt *vertMask, const PetscInt currentLevel);
+      PetscInt FFM_VertexBased_SolveModified(const PetscInt dim, const PetscInt minVerts, const PetscReal x0[], const PetscInt nVert, PetscInt verts[], PetscScalar *updatedVertex, PetscScalar *lsArray, PetscReal *updatedLS, const PetscInt *vertMask, const PetscInt currentLevel);
+	  PetscInt FFM_VertexBased_ModifiedGreenGauss(const PetscInt dim, PetscInt id_potential, PetscInt id_base, PetscScalar *updatedVertex, PetscScalar *xCoord, PetscScalar *yCoord, const PetscInt *vertMask, const PetscInt currentLevel, PetscScalar *lsArray, PetscReal *updatedLS);
+	  PetscInt FFM_VertexBased_RBF(const PetscInt dim, PetscInt id_potential, PetscInt id_base, PetscScalar *updatedVertex, PetscScalar *xCoord, PetscScalar *yCoord, const PetscInt *vertMask, const PetscInt currentLevel, PetscScalar *lsArray, PetscReal *updatedLS);
+	  PetscInt FFM_VertexBased_GMGG(const PetscInt dim, PetscInt id_potential, PetscInt id_base, PetscScalar *updatedVertex, PetscScalar *xCoord, PetscScalar *yCoord, const PetscInt *vertMask, const PetscInt currentLevel, PetscScalar *lsArray, PetscReal *updatedLS);
 
       std::shared_ptr<ablate::levelSet::GaussianConvolution> convolution = nullptr;
 
@@ -130,7 +134,32 @@ namespace ablate::levelSet {
 
       // Given a cell-centered VOF field compute the level-set field
       void ToLevelSet(DM vofDM, Vec vofVec, const ablate::domain::Field vofField);
-
+
+      PetscScalar LScircle(const PetscReal* x, const PetscInt dim);
+      PetscBool CutCellfromLS(DM aux_dm, const PetscInt point, const ablate::domain::Field *levelSetField, Vec auxVector);
+      void arbit_interface(DM dm, const ablate::domain::Field levelSetField, Vec auxVec);
+      PetscScalar EllipseEq(const PetscReal* x, PetscScalar theta, PetscScalar a, PetscScalar b, PetscScalar* df, PetscScalar* d2f);
+      PetscScalar LSellipse(const PetscReal* x, const PetscInt dim);
+      PetscScalar Newton(const PetscReal* x, void* ctx);
+      PetscScalar CassiniOvalEq(const PetscReal* x, PetscScalar theta, PetscScalar a, PetscScalar c, PetscInt branch, PetscScalar* df, PetscScalar* d2f);
+	  PetscScalar LScassinioval(const PetscReal* x, const PetscInt dim);
+	  PetscScalar NewtonCassiniOval(const PetscReal* x, void* ctx);
+	  PetscScalar StarEq(const PetscReal* x, PetscScalar theta, PetscScalar a, PetscScalar b, PetscScalar* df, PetscScalar* d2f);
+	  PetscScalar LSstar(const PetscReal* x, const PetscInt dim);
+	  PetscScalar NewtonStar(const PetscReal* x, void* ctx, PetscScalar* theta_star);
+	  PetscScalar LStwoCircles(const PetscReal* x, const PetscInt dim);
+  };
+
+  struct EllipseCtx {
+    PetscScalar a, b;
+  };
+
+  struct CassiniOvalCtx {
+	  PetscScalar a, e, c;
+  };
+
+  struct StarCtx {
+	  PetscScalar a, b;
   };
 
 

src/levelSet/levelSetUtilities.cpp

@@ -1076,7 +1076,7 @@ void ablate::levelSet::Utilities::Reinitialize(
 
 #ifdef saveData
   char fname[255];
-  sprintf(fname, "vof0_%03ld.txt", saveIter);
+  sprintf(fname, "vof0_%03d.txt", saveIter);
   SaveCellData(solDM, solVec, fname, vofField, 1, subDomain);
 #endif
 
@@ -1131,7 +1131,7 @@ VecGetArray(workVec, &workArray);
 
 
 #ifdef saveData
-  sprintf(fname, "vof1_%03ld.txt", saveIter);
+  sprintf(fname, "vof1_%03d.txt", saveIter);
   SaveCellData(auxDM, workVec, fname, vofField, 1, subDomain);
 #endif
 
@@ -1154,14 +1154,14 @@ VecGetArray(workVec, &workArray);
 
 #ifdef saveData
 {
-  sprintf(fname, "mask0_%03ld.txt", saveIter);
+  sprintf(fname, "mask0_%03d.txt", saveIter);
   FILE *f1 = fopen(fname, "w");
   for (PetscInt c = cellRange.start; c < cellRange.end; ++c) {
     const PetscInt cell = cellRange.GetPoint(c);
     PetscReal x[dim];
     DMPlexComputeCellGeometryFVM(solDM, cell, NULL, x, NULL) >> ablate::utilities::PetscUtilities::checkError;
     for (PetscInt d = 0; d < dim; ++d) fprintf(f1, "%+f\t", x[d]);
-    fprintf(f1, "%ld\n", cellMask[c]);
+    fprintf(f1, "%d\n", cellMask[c]);
   }
   fclose(f1);
 }
@@ -1268,17 +1268,17 @@ for (PetscInt c = cellRange.start; c < cellRange.end; ++c) {
 
 #ifdef saveData
 {
-  sprintf(fname, "mask1_%03ld.txt", saveIter);
+  sprintf(fname, "mask1_%03d.txt", saveIter);
   FILE *f1 = fopen(fname, "w");
   for (PetscInt c = cellRange.start; c < cellRange.end; ++c) {
     const PetscInt cell = cellRange.GetPoint(c);
     PetscReal x[dim];
     DMPlexComputeCellGeometryFVM(solDM, cell, NULL, x, NULL) >> ablate::utilities::PetscUtilities::checkError;
     for (PetscInt d = 0; d < dim; ++d) fprintf(f1, "%+f\t", x[d]);
-    fprintf(f1, "%ld\n", cellMask[c]);
+    fprintf(f1, "%d\n", cellMask[c]);
   }
   fclose(f1);
-  sprintf(fname, "cellNormal0_%03ld.txt", saveIter);
+  sprintf(fname, "cellNormal0_%03d.txt", saveIter);
   SaveCellData(auxDM, auxVec, fname, cellNormalField, dim, subDomain);
 }
 #endif
@@ -1361,7 +1361,7 @@ for (PetscInt c = cellRange.start; c < cellRange.end; ++c) {
 
 
 #ifdef saveData
-  sprintf(fname, "ls1_%03ld.txt", saveIter);
+  sprintf(fname, "ls1_%03d.txt", saveIter);
   SaveVertexData(auxDM, auxVec, fname, lsField, 1, subDomain);
 #endif
 
@@ -1498,17 +1498,17 @@ for (PetscInt c = cellRange.start; c < cellRange.end; ++c) {
 
 #ifdef saveData
 {
-  sprintf(fname, "mask2_%03ld.txt", saveIter);
+  sprintf(fname, "mask2_%03d.txt", saveIter);
   FILE *f1 = fopen(fname, "w");
   for (PetscInt c = cellRange.start; c < cellRange.end; ++c) {
     const PetscInt cell = cellRange.GetPoint(c);
     PetscReal x[dim];
     DMPlexComputeCellGeometryFVM(solDM, cell, NULL, x, NULL) >> ablate::utilities::PetscUtilities::checkError;
     for (PetscInt d = 0; d < dim; ++d) fprintf(f1, "%+f\t", x[d]);
-    fprintf(f1, "%ld\n", cellMask[c]);
+    fprintf(f1, "%d\n", cellMask[c]);
   }
   fclose(f1);
-  sprintf(fname, "ls2_%03ld.txt", saveIter);
+  sprintf(fname, "ls2_%03d.txt", saveIter);
   SaveVertexData(auxDM, auxVec, fname, lsField, 1, subDomain);
 }
 #endif
@@ -1625,7 +1625,7 @@ for (PetscInt c = cellRange.start; c < cellRange.end; ++c) {
   }
 
 #ifdef saveData
-  sprintf(fname, "ls3_%03ld.txt", saveIter);
+  sprintf(fname, "ls3_%03d.txt", saveIter);
   SaveVertexData(auxDM, auxVec, fname, lsField, 1, subDomain);
 #endif
 
@@ -1641,7 +1641,7 @@ for (PetscInt c = cellRange.start; c < cellRange.end; ++c) {
   }
 
 #ifdef saveData
-  sprintf(fname, "mask3_%03ld.txt", saveIter);
+  sprintf(fname, "mask3_%03d.txt", saveIter);
   SaveCellData(auxDM, workVec, fname, vofField, 1, subDomain);
 #endif
 
@@ -1664,7 +1664,7 @@ for (PetscInt c = cellRange.start; c < cellRange.end; ++c) {
 
   subDomain->UpdateAuxLocalVector();
 #ifdef saveData
-  sprintf(fname, "curv0_%03ld.txt", saveIter);
+  sprintf(fname, "curv0_%03d.txt", saveIter);
   SaveCellData(auxDM, auxVec, fname, curvID, 1, subDomain);
 #endif
 
@@ -1711,7 +1711,7 @@ for (PetscInt c = cellRange.start; c < cellRange.end; ++c) {
 
 
 #ifdef saveData
-  sprintf(fname, "vertH0_%03ld.txt", saveIter);
+  sprintf(fname, "vertH0_%03d.txt", saveIter);
   SaveVertexData(auxDM, workVec, fname, lsField, 1, subDomain);
 #endif
 
@@ -1800,7 +1800,7 @@ for (PetscInt c = cellRange.start; c < cellRange.end; ++c) {
 
 
 #ifdef saveData
-  sprintf(fname, "vertH1_%03ld.txt", saveIter);
+  sprintf(fname, "vertH1_%03d.txt", saveIter);
   SaveVertexData(auxDM, workVec, fname, lsField, 1, subDomain);
 #endif
 
@@ -1853,7 +1853,7 @@ for (PetscInt c = cellRange.start; c < cellRange.end; ++c) {
 
 
 #ifdef saveData
-  sprintf(fname, "vertH2_%03ld.txt", saveIter);
+  sprintf(fname, "vertH2_%03d.txt", saveIter);
   SaveVertexData(auxDM, workVec, fname, lsField, 1, subDomain);
 #endif
 
@@ -1887,9 +1887,9 @@ for (PetscInt c = cellRange.start; c < cellRange.end; ++c) {
   subDomain->UpdateAuxLocalVector();
 
 #ifdef saveData
-  sprintf(fname, "cellH0_%03ld.txt", saveIter);
+  sprintf(fname, "cellH0_%03d.txt", saveIter);
   SaveVertexData(auxDM, workVec, fname, lsField, 1, subDomain);
-  sprintf(fname, "cellNormal1_%03ld.txt", saveIter);
+  sprintf(fname, "cellNormal1_%03d.txt", saveIter);
   SaveCellData(auxDM, auxVec, fname, cellNormalField, dim, subDomain);
 
 #endif

src/utilities/petscSupport.cpp

@@ -109,7 +109,7 @@
       PetscReal x[dim];
       PetscCall(DMPlexComputeCellGeometryFVM(dm, cell, NULL, x, NULL));
       printf("plot(%f,%f,'r*'); %% Cell\n", x[0], x[1]);
-printf("%ld\n", sharedFace);
+printf("%d\n", sharedFace);
 //      PetscCall(DMPlexComputeCellGeometryFVM(dm, sharedFace, NULL, x, NULL));
 //      printf("plot(%f,%f,'r*'); %% Shared Face\n", x[0], x[1]);
       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "DMPlexGetForwardCell detected that a face is shared between %" PetscInt_FMT" cells.", nPoints);

src/utilities/petscSupport.hpp

@@ -51,7 +51,7 @@ PetscErrorCode DMPlexFaceCentroidOutwardAreaNormal(DM dm, PetscInt cell, PetscIn
  * maxLevels - Number of neighboring cells/vertices to check
  * maxDist - Maximum distance to include
  * numberCells - The number of cells/vertices to return.
- * useSharedFace - Return cells/vertices which share a common face (PETSC_TRUE) or a shared vertex (PETSC_FALSE)
+ * useSharedFace - Return cells/vertices which share a common face (PETSC_TRUE) or a shared edge (PETSC_FALSE) 
  * returnVertices - Return vertices surrounding the center cell (PETSC_TRUE) or cells surrounding the center cell (PETSC_FALSE)
  * nCells - Number of neighboring cells/vertices
  * cells - The list of neighboring cell/vertices IDs

tests/regressionTests/inputs/exampleRegressionTest/_exampleIntegrationTest/exampleRegressionTest.yaml

@@ -0,0 +1,18 @@
+---
+test:
+  name: exampleRegressionTest
+  assert: inputs/exampleRegressionTest/expectedOutput.txt
+environment:
+  title: _exampleIntegrationTest
+  tagDirectory: false
+arguments: {}
+timestepper:
+  arguments: {}
+  domain: !<!ablate::domain::BoxMesh>
+    name: simpleBoxField
+    faces: [12, 12]
+    lower: [0, 0]
+    upper: [1, 1]
+    modifiers: []
+    fields: []
+solvers: []