Added GenerateMesh.h in basis folder

dftfeDevelopers · Jun 15, 2024 · 6a4cb05 · 6a4cb05
1 parent 114b143
commit 6a4cb05
Show file tree

Hide file tree

Showing 29 changed files with 1,349 additions and 32 deletions.
diff --git a/install.py b/install.py
@@ -2,10 +2,10 @@
 import os
 import sys
 import textwrap
-opts_dict={'build_dir' : './build'}
+opts_dict={'build_dir' : './build', 'n' : ''}
 
 def getUsageMsg():
-    msg = '''The correct usage is python install.py '''\
+    msg = '''The correct usage is python install.py [--n=<numprocs>]\n'''\
           '''[--build_dir=/path/to/build/directory]\n'''\
           '''[--src_dir=/path/to/source/directory]\n'''\
           '''The optional [--build_dir=/path/to/build/directory] specifies '''\
@@ -79,6 +79,7 @@ def updateOptsDictFromCommandLine(strings):
             y = y + " " + x
             print(x, end=" ")
 
+        nprocs = opts_dict['n']
         print()
         os.system(y)
-        os.system('make -j')
+        os.system('make -j '+nprocs)
diff --git a/src/basis/CMakeLists.txt b/src/basis/CMakeLists.txt
@@ -31,6 +31,7 @@ set(DFT-EFE-BASIS-SOURCES
 	EFEBasisDataStorageDealii.cpp
 	EFEBasisDofHandlerDealii.cpp
 	EFEConstraintsLocalDealii.cpp
+	GenerateMesh.cpp
 	)
 
 add_library(dft-efe-basis SHARED ${DFT-EFE-BASIS-SOURCES})

diff --git a/src/basis/Defaults.cpp b/src/basis/Defaults.cpp
@@ -38,6 +38,6 @@ namespace dftefe
     const double    L2ProjectionDefaults::ABSOLUTE_TOL           = 1e-13;
     const double    L2ProjectionDefaults::RELATIVE_TOL           = 1e-14;
     const double    L2ProjectionDefaults::DIVERGENCE_TOL         = 1e6;
-
+    const size_type GenerateMeshDefaults::MAX_REFINEMENT_STEPS   = 40;
   } // end of namespace basis
 } // end of namespace dftefe
diff --git a/src/basis/Defaults.h b/src/basis/Defaults.h
@@ -69,6 +69,17 @@ namespace dftefe
 
     }; // end of class L2ProjectionDefaults
 
+    class GenerateMeshDefaults
+    {
+    public:
+      //
+      // The maximum refinement steps used for adaptive mesh generation
+      // algorithm.
+      //
+      static const size_type MAX_REFINEMENT_STEPS;
+
+    }; // end of class L2ProjectionDefaults
+
   } // end of namespace basis
 } // end of namespace dftefe
 #endif // dftefebasisDefaults_h
diff --git a/src/basis/EnrichmentClassicalInterfaceSpherical.t.cpp b/src/basis/EnrichmentClassicalInterfaceSpherical.t.cpp
@@ -147,7 +147,9 @@ namespace dftefe
           fieldName,
           minbound,
           maxbound,
-          d_triangulation->maxCellDiameter(),
+          d_triangulation->maxElementLength(),
+          d_triangulation->getDomainVectors(),
+          d_triangulation->getPeriodicFlags(),
           cellVerticesVector,
           comm);
 
@@ -465,6 +467,8 @@ namespace dftefe
           minbound,
           maxbound,
           0,
+          triangulation->getDomainVectors(),
+          d_triangulation->getPeriodicFlags(),
           cellVerticesVector,
           comm);
 

diff --git a/src/basis/EnrichmentIdsPartition.h b/src/basis/EnrichmentIdsPartition.h
@@ -84,6 +84,8 @@ namespace dftefe
         const std::vector<double> &                  minbound,
         const std::vector<double> &                  maxbound,
         double                                       additionalCutoff,
+        const std::vector<utils::Point> &            globalDomainBoundVec,
+        const std::vector<bool> &                    isPeriodicFlags,
         const std::vector<std::vector<utils::Point>> &cellVerticesVector,
         const utils::mpi::MPIComm &                   comm);
 

diff --git a/src/basis/EnrichmentIdsPartition.t.cpp b/src/basis/EnrichmentIdsPartition.t.cpp
@@ -195,8 +195,33 @@ namespace dftefe
         const std::vector<double> &                  minbound,
         const std::vector<double> &                  maxbound,
         double                                       additionalCutoff,
-        const std::vector<std::vector<utils::Point>> &cellVerticesVector)
+        const std::vector<bool> &                    isPeriodicFlags,
+        const std::vector<std::vector<utils::Point>> &cellVerticesVector,
+        const utils::mpi::MPIComm &                   comm)
       {
+        std::vector<double> minboundGlobalDomain(dim, 0.),
+          maxboundGlobalDomain(dim, 0.);
+        if (!(std::all_of(isPeriodicFlags.begin(),
+                          isPeriodicFlags.end(),
+                          [](bool v) { return v; })))
+          {
+            int err = utils::mpi::MPIAllreduce<utils::MemorySpace::HOST>(
+              minbound.data(),
+              minboundGlobalDomain.data(),
+              minbound.size(),
+              utils::mpi::MPIDouble,
+              utils::mpi::MPIMin,
+              comm);
+
+            err = utils::mpi::MPIAllreduce<utils::MemorySpace::HOST>(
+              maxbound.data(),
+              maxboundGlobalDomain.data(),
+              maxbound.size(),
+              utils::mpi::MPIDouble,
+              utils::mpi::MPIMax,
+              comm);
+          }
+
         std::vector<size_type> newAtomIds = atomIdsPartition->newAtomIds();
         std::vector<double>    minCellBound;
         std::vector<double>    maxCellBound;
@@ -264,6 +289,28 @@ namespace dftefe
 
                         DFTEFE_Assert(b >= a);
                         DFTEFE_Assert(d >= c);
+
+                        // check that enrichment functions do not spill bounding
+                        // box for non-periodic cases.
+                        if (!isPeriodicFlags[k] &&
+                            !(c - additionalCutoff > minboundGlobalDomain[k] &&
+                              d + additionalCutoff < maxboundGlobalDomain[k]))
+                          {
+                            std::string msg =
+                              "The enrichment functions may spill to a"
+                              " non-periodic face of the triangulation domain which is not allowed."
+                              " Increase the "
+                              " domain boundary or reduce to the ball radius of the enrichment "
+                              " function cutoff. ";
+                            if (additionalCutoff != 0)
+                              msg +=
+                                "Recommended domain boundary increase, if wanted, can be by " +
+                                std::to_string(additionalCutoff) +
+                                " bohr on each side of origin.";
+                            utils::throwException<utils::InvalidArgument>(false,
+                                                                          msg);
+                          }
+
                         if (!((c < a && d < a) || (c > b && d > b)))
                           flag = true;
                         else
@@ -366,9 +413,34 @@ namespace dftefe
       const std::vector<double> &                   minbound,
       const std::vector<double> &                   maxbound,
       double                                        additionalCutoff,
+      const std::vector<utils::Point> &             globalDomainBoundVec,
+      const std::vector<bool> &                     isPeriodicFlags,
       const std::vector<std::vector<utils::Point>> &cellVerticesVector,
       const utils::mpi::MPIComm &                   comm)
     {
+      utils::throwException<utils::InvalidArgument>(
+        !(std::any_of(isPeriodicFlags.begin(),
+                      isPeriodicFlags.end(),
+                      [](bool v) { return v; })),
+        "EnrichmentIdsPartition can only handle non-periodic boundary conditions."
+        " Contact Developers to get it extended to periodic systems.");
+
+      double sum = 0.0;
+      for (size_type i = 0; i < globalDomainBoundVec.size(); i++)
+        {
+          for (size_type j = 0; j < dim; j++)
+            {
+              if (i != j)
+                sum += globalDomainBoundVec[i][j];
+            }
+        }
+      utils::throwException<utils::InvalidArgument>(
+        sum < 1e-12,
+        "EnrichmentIdsPartition can only handle orthogonal domains with cartesian"
+        "coordinate domain vectors {(1,0,0), (0,1,0, (0,0,1)}. Contact Developers"
+        " to get it extended to non orthogonal systems.");
+      // Note this class cannot handle rotated orthogonal domain also...
+
       std::vector<double>    rCutoffMax;
       std::vector<size_type> atomIds;
       rCutoffMax.resize(atomSymbol.size(), 0.);
@@ -396,12 +468,15 @@ namespace dftefe
         atomSymbol,
         fieldName,
         comm);
+
       EnrichmentIdsPartitionInternal::getLocalEnrichmentIds<dim>(
         d_locallyOwnedEnrichmentIds,
         d_newAtomIdToEnrichmentIdOffset,
         atomIdsPartition);
+
       EnrichmentIdsPartitionInternal::getOverlappingAtomIdsInBox<dim>(
         atomIds, rCutoffMax, atomCoordinates, minbound, maxbound);
+
       EnrichmentIdsPartitionInternal::getOverlappingEnrichmentIdsInCells<dim>(
         d_overlappingEnrichmentIdsInCells,
         atomIds,
@@ -414,7 +489,10 @@ namespace dftefe
         minbound,
         maxbound,
         additionalCutoff,
-        cellVerticesVector);
+        isPeriodicFlags,
+        cellVerticesVector,
+        comm);
+
       EnrichmentIdsPartitionInternal::getGhostEnrichmentIds<dim>(
         atomIdsPartition,
         d_enrichmentIdsInProcessor,

diff --git a/src/basis/FECellBase.h b/src/basis/FECellBase.h
@@ -55,6 +55,12 @@ namespace dftefe
       virtual void
       clearRefineFlag() = 0;
 
+      virtual double
+      minimumVertexDistance() const = 0;
+
+      virtual double
+      distanceToUnitCell(dftefe::utils::Point &parametricPoint) const = 0;
+
       virtual void
       setCoarsenFlag() = 0;
 

diff --git a/src/basis/FECellDealii.h b/src/basis/FECellDealii.h
@@ -76,6 +76,12 @@ namespace dftefe
       void
       clearRefineFlag() override;
 
+      double
+      minimumVertexDistance() const override;
+
+      double
+      distanceToUnitCell(dftefe::utils::Point &parametricPoint) const override;
+
       void
       setCoarsenFlag() override;
 

diff --git a/src/basis/FECellDealii.t.cpp b/src/basis/FECellDealii.t.cpp
@@ -137,6 +137,23 @@ namespace dftefe
       d_dealiiFECellIter->clear_refine_flag();
     }
 
+    template <unsigned int dim>
+    double
+    FECellDealii<dim>::minimumVertexDistance() const
+    {
+      return d_dealiiFECellIter->minimum_vertex_distance();
+    }
+
+    template <unsigned int dim>
+    double
+    FECellDealii<dim>::distanceToUnitCell(
+      dftefe::utils::Point &parametricPoint) const
+    {
+      dealii::Point<dim, double> dealiiPoint;
+      convertToDealiiPoint<dim>(parametricPoint, dealiiPoint);
+      return dealii::GeometryInfo<dim>::distance_to_unit_cell(dealiiPoint);
+    }
+
     template <unsigned int dim>
     void
     FECellDealii<dim>::setCoarsenFlag()