Mesh volume conservation

hoomd/MeshDefinition.cc

@@ -28,10 +28,15 @@ MeshDefinition::MeshDefinition(std::shared_ptr<SystemDefinition> sysdef, unsigne
           std::shared_ptr<TriangleData>(new TriangleData(m_sysdef->getParticleData(), n_types)))
 
     {
+    // allocate the max number of neighbors per type allowed
+    GlobalArray<unsigned int> globalN(n_types, m_sysdef->getParticleData()->getExecConf());
+    m_globalN.swap(globalN);
+    TAG_ALLOCATION(m_globalN);
     }
 
 void MeshDefinition::setTypes(pybind11::list types)
     {
+    m_globalN.resize(len(types));
     for (unsigned int i = 0; i < len(types); i++)
         {
         m_meshbond_data->setTypeName(i, types[i].cast<string>());
@@ -106,13 +111,20 @@ void MeshDefinition::setTriangulationData(pybind11::dict triangulation)
     triangle_data.resize(static_cast<unsigned int>(len_triang));
     TriangleData::members_t triangle_new;
 
+    ArrayHandle<unsigned int> h_globalN(m_globalN, access_location::host, access_mode::overwrite);
+
+    for (unsigned int i = 0; i < m_meshtriangle_data->getNTypes(); i++)
+        h_globalN.data[i] = 0;
+
     for (size_t i = 0; i < len_triang; i++)
         {
         triangle_new.tag[0] = ptr1[i * 3];
         triangle_new.tag[1] = ptr1[i * 3 + 1];
         triangle_new.tag[2] = ptr1[i * 3 + 2];
         triangle_data.groups[i] = triangle_new;
         triangle_data.type_id[i] = ptr2[i];
+
+        h_globalN.data[triangle_data.type_id[i]] += 1;
         }
 
     m_meshtriangle_data = std::shared_ptr<TriangleData>(

hoomd/MeshDefinition.h

@@ -9,6 +9,8 @@
 #error This header cannot be compiled by nvcc
 #endif
 
+#include "hoomd/GlobalArray.h"
+
 #include "BondedGroupData.h"
 #include "MeshGroupData.h"
 #include "SystemDefinition.h"
@@ -85,6 +87,11 @@ class PYBIND11_EXPORT MeshDefinition
         return triangles.getSize();
         }
 
+    const GlobalArray<unsigned int>& getPerTypeSize() const
+        {
+        return m_globalN;
+        }
+
     void setTypes(pybind11::list types);
 
     BondData::Snapshot getBondData();
@@ -96,6 +103,7 @@ class PYBIND11_EXPORT MeshDefinition
     void setTriangulationData(pybind11::dict triangulation);
 
     private:
+    GlobalArray<unsigned int> m_globalN;
     std::shared_ptr<SystemDefinition>
         m_sysdef; //!< System definition later needed for dynamic bonding
     std::shared_ptr<MeshBondData> m_meshbond_data;     //!< Bond data for the mesh

hoomd/md/CMakeLists.txt

@@ -54,6 +54,7 @@ set(_md_sources module-md.cc
                    Thermostat.cc
                    TwoStepNVTAlchemy.cc
                    WallData.cc
+		   VolumeConservationMeshForceCompute.cc
                    ZeroMomentumUpdater.cc
                    )
 
@@ -206,6 +207,9 @@ set(_md_headers ActiveForceComputeGPU.h
                 TwoStepConstantPressure.h
                 AlchemostatTwoStep.h
                 TwoStepNVTAlchemy.h
+		VolumeConservationMeshForceCompute.h
+		VolumeConservationMeshForceComputeGPU.h
+                AlchemostatTwoStep.h
                 WallData.h
                 ZeroMomentumUpdater.h
                 )
@@ -241,6 +245,7 @@ list(APPEND _md_sources ActiveForceComputeGPU.cc
                            TableDihedralForceComputeGPU.cc
                            TwoStepBDGPU.cc
                            TwoStepLangevinGPU.cc
+			   VolumeConservationMeshForceComputeGPU.cc
                            TwoStepConstantVolumeGPU.cc
                            TwoStepConstantPressureGPU.cc
                            MuellerPlatheFlowGPU.cc
@@ -283,6 +288,7 @@ set(_md_cu_sources ActiveForceComputeGPU.cu
                       TwoStepConstantVolumeGPU.cu
                       TwoStepNVEGPU.cu
                       TwoStepRATTLENVEGPU.cu
+		      VolumeConservationMeshForceComputeGPU.cu
                       MuellerPlatheFlowGPU.cu
                       CosineSqAngleForceGPU.cu
                       )