simplify the assembly process

src/mumfim/analysis/ElementalSystem.cc

@@ -1,38 +1,66 @@
 #include "ElementalSystem.h"
+#include <apfNumbering.h>
+#include <apfField.h>
+
 namespace amsi
 {
-  ElementalSystem::ElementalSystem(apf::Field * field, int o)
-    : apf::Integrator(o)
-    , Ke()
-    , fe()
-    , me(NULL)
-    , e(NULL)
-    , f(field)
-    , nedofs(0)
-    , nenodes(0)
-    , num_field_components(0)
+  // this function works for both scalar and vector fields
+  static void GetNodalFieldValuesAndNumbers(apf::Field* field,
+                                            apf::Numbering* numbering,
+                                            apf::MeshEntity* mesh_entity, std::vector<int>& dof_numbers, std::vector<double> & values) {
+    auto nnodes = field->countNodesOn(mesh_entity);
+    const auto ncomponents = apf::countComponents(field);
+    values.resize(ncomponents*nnodes);
+    dof_numbers.resize(ncomponents*nnodes);
+    for(int node=0; node<nnodes; ++node) {
+      apf::getComponents(field, mesh_entity, node, &values[node*ncomponents]);
+      for(int component=0; component<ncomponents; ++component) {
+        dof_numbers[node*ncomponents+component] = apf::getNumber(numbering,mesh_entity,node, component);
+      }
+    }
+  }
+  ElementalSystem::ElementalSystem(apf::Field * field, apf::Numbering* numbering, int o)
+      : apf::Integrator(o)
+      , Ke()
+      , fe()
+      , me(NULL)
+      , e(NULL)
+      , f(field)
+      , numbering_(numbering)
+      , nedofs(0)
+      , nenodes(0)
+      , num_field_components(0)
   {
     num_field_components = apf::countComponents(f);
   }
-  // if the new mesh elemenent is on a mesh entity with the same type as the previous mesh element was, we don't need to reallocate the memory for the element matrix and vector (just zero it), otherwise we need to resize it
+
+  // if the new mesh elemenent is on a mesh entity with the same type as the
+  // previous mesh element was, we don't need to reallocate the memory for the
+  // element matrix and vector (just zero it), otherwise we need to resize it
   void ElementalSystem::inElement(apf::MeshElement * ME)
   {
-    e = apf::createElement(f,me);
+    me = ME;
+    e = apf::createElement(f, me);
     nenodes = apf::countNodes(e);
     int new_nedofs = nenodes * num_field_components;
     bool reallocate = nedofs != new_nedofs;
     nedofs = new_nedofs;
-    if(reallocate)
+    GetNodalFieldValuesAndNumbers(f, numbering_, getMeshEntity(ME),
+                                  field_numbers_, field_values_);
+        // TODO: remove this check. login in DynamicArray will account for this
+    // Note: the reallocation routine in DynamicArray will cause a lot of reallocations
+    if (reallocate)
     {
-      Ke.setSize(nedofs,nedofs);
+      Ke.setSize(nedofs, nedofs);
       fe.setSize(nedofs);
     }
     Ke.zero();
     fe.zero();
   }
+
   void ElementalSystem::outElement()
   {
     apf::destroyElement(e);
     e = nullptr;
   }
-}
+}  // namespace amsi

src/mumfim/analysis/ElementalSystem.h

@@ -2,6 +2,7 @@
 #define ELEMENTALSYSTEM_H_
 #include <apf.h>
 #include <apfDynamicMatrix.h>
+#include <apfNumbering.h>
 namespace amsi
 {
   class ElementalSystem : public apf::Integrator
@@ -12,21 +13,25 @@ namespace amsi
     apf::MeshElement * me;
     apf::Element * e;
     apf::Field * f;
+    apf::Numbering* numbering_;
     int nedofs;
     int nenodes;
     int num_field_components;
+    std::vector<double> field_values_;
+    std::vector<int> field_numbers_;
   public:
-    ElementalSystem(apf::Field * f, int o);
+    ElementalSystem(apf::Field * f, apf::Numbering* numbering, int o);
     virtual void inElement(apf::MeshElement *);
     virtual void outElement();
     virtual void parallelReduce() {};
     virtual bool includesBodyForces() {return false;}
     virtual int numElementalDOFs() {return nedofs;}
-    // This shouldn't be used because it causes segfaults when
-    // used improperly (which is almost always!
-    // virtual apf::Element * getElement() {return e;}
     virtual apf::DynamicMatrix& getKe() {return Ke;}
     virtual apf::DynamicVector& getfe() {return fe;}
+
+    // TODO have getFieldValues and getFieldNumbers return mdspan
+    [[nodiscard]] const std::vector<double>& getFieldValues() { return field_values_; }
+    [[ nodiscard ]] const std::vector<int> & getFieldNumbers() { return field_numbers_; }
     virtual apf::Field * getField() {return f;}
   };
 }

src/mumfim/analysis/amsiFEA.cc

@@ -1,5 +1,6 @@
 #include "amsiFEA.h"
 
+#include <apfFieldData.h>
 #include <apfNumbering.h>
 #include <apfShape.h>
 #include <maShape.h>
@@ -50,8 +51,9 @@ namespace amsi
       {
         if (apf_mesh->isOwned(mesh_ent))
         {
-          apf::FieldShape * fs = apf::getShape(apf_primary_field);
-          int num_nodes = fs->countNodesOn(apf_mesh->getType(mesh_ent));
+          apf::FieldShape * shape_function = apf::getShape(apf_primary_field);
+          int num_nodes =
+              shape_function->countNodesOn(apf_mesh->getType(mesh_ent));
           for (int jj = 0; jj < num_nodes; jj++)
           {
             apf::Vector3 disp;
@@ -128,17 +130,18 @@ namespace amsi
 
   static ModelDefinition GetModelDefinitionComponent(
       const mt::CategoryNode & analysis_case,
-      const std::string & analysis_name, const std::string& component_name)
+      const std::string & analysis_name,
+      const std::string & component_name)
   {
-    auto* component = mt::GetCategoryByType(&analysis_case, component_name);
-    auto* temp = mt::GetCategoryByType(component, analysis_name);
+    auto * component = mt::GetCategoryByType(&analysis_case, component_name);
+    auto * temp = mt::GetCategoryByType(component, analysis_name);
     component = (temp == nullptr) ? component : temp;
-    if(component == nullptr) {
+    if (component == nullptr)
+    {
       throw mumfim::mumfim_error(component_name + " not defined");
     }
-    return {
-        .associated = mt::AssociateModel(component),
-        .unassociated = std::make_shared<mt::CategoryNode>(*component)};
+    return {.associated = mt::AssociateModel(component),
+            .unassociated = std::make_shared<mt::CategoryNode>(*component)};
   }
 
   FEAStep::FEAStep(apf::Mesh * mesh,
@@ -187,114 +190,61 @@ namespace amsi
     offset = first_local_dof;
   }
 
-  template <typename NODE_TYPE>
-  static void AssembleDOFs(LAS* las, int num_elemental_dofs, int* dof_numbers,
-                             const NODE_TYPE* node_values, double* Ke, double* fe,
-                             bool includes_body_forces, int analysis_dim)
+  static void AssembleDOFs(LAS * las,
+                           const std::vector<int> & dof_numbers,
+                           const std::vector<double> & dof_values,
+                           apf::DynamicMatrix & Ke,
+                           apf::DynamicVector & fe,
+                           bool includes_body_forces,
+                           int ncomponents)
   {
-    if (Ke != NULL) {
-      if (!includes_body_forces) {
-        double* bf = new double[num_elemental_dofs]();
-        for (int ii = 0; ii < num_elemental_dofs; ii++) {
-          const int& global_i = dof_numbers[ii];
-          // this is the isFixed function from apfFunctions
-          // which is different from the sim query is fixed function!
-          if (!isFixed(global_i)) {
-            for (int jj = 0; jj < num_elemental_dofs; jj++) {
-              const double& val = Ke[ii * num_elemental_dofs + jj];
-              double j_val = node_values[jj / analysis_dim][jj % analysis_dim];
-              if (j_val != 0.0) bf[ii] += -val * j_val;
-            }
-          }
-        }
-        las->AddToVector(num_elemental_dofs, dof_numbers, &bf[0]);
-        delete[] bf;
-      }
-      las->AddToMatrix(num_elemental_dofs, dof_numbers, num_elemental_dofs,
-                       dof_numbers, &Ke[0]);
-    }
-    /// Modification of fe to correctly account for nonzero dirichlet boundary
-    /// conditions
-    double* dirichletValue = new double[num_elemental_dofs]();
-    for (int ii = 0; ii < num_elemental_dofs; ii++) {
-      if (dof_numbers[ii] < 0)
-        dirichletValue[ii] = node_values[ii / analysis_dim][ii % analysis_dim];
+    if (!includes_body_forces)
+    {
+      throw mumfim::mumfim_error(
+          "Integrator that does not bring int it's own body forces is no "
+          "longer supported");
     }
-    double* dfe = new double[num_elemental_dofs]();
+    const auto num_elemental_dofs = fe.size();
+    las->AddToMatrix(num_elemental_dofs, dof_numbers.data(), num_elemental_dofs,
+                     dof_numbers.data(), &Ke(0, 0));
     for (int ii = 0; ii < num_elemental_dofs; ii++)
+    {
+      // if the dof_number is negative than that particular dof is fixed
+      // so set the value to zero if it's not fixed and to the real value if
+      // it is fixed.
+      double dirichlet_value = dof_numbers[ii] < 0 ? dof_values[ii] : 0.0;
       for (int jj = 0; jj < num_elemental_dofs; jj++)
-        dfe[ii] =
-            dfe[ii] + Ke[ii * num_elemental_dofs + jj] * dirichletValue[ii];
-    for (int ii = 0; ii < num_elemental_dofs; ii++)
-      fe[ii] = fe[ii] + dfe[ii];
-    delete[] dirichletValue;
-    delete[] dfe;
-    las->AddToVector(num_elemental_dofs, dof_numbers, &fe[0]);
-  }
-
-  template <>
-  void AssembleDOFs<double>(LAS* las, int num_elemental_dofs, int* dof_numbers,
-                                            const double* node_values, double* Ke, double* fe,
-                                            bool includes_body_forces, int analysis_dim)
-  {
-    if (Ke == NULL) {
-      throw mumfim::mumfim_error("Ke is null going into AssembleDOFs");
-    }
-    if (!includes_body_forces) {
-      throw mumfim::mumfim_error("Scalar does not handle case with !includes_body_force");
-    }
-
-    las->AddToMatrix(num_elemental_dofs, dof_numbers, num_elemental_dofs,
-                     dof_numbers, &Ke[0]);
-
-    /// Modification of fe to correctly account for nonzero dirichlet boundary
-    /// conditions
-    std::vector<double> dirichletValue(num_elemental_dofs, 0.0);
-    for (int ii = 0; ii < num_elemental_dofs; ii++) {
-      if (dof_numbers[ii] < 0) dirichletValue[ii] = node_values[ii];
+      {
+        fe[ii] += Ke(ii, jj) * dirichlet_value;
+      }
     }
-    std::vector<double> dfe(num_elemental_dofs, 0.0);
-    for (int ii = 0; ii < num_elemental_dofs; ii++)
-      for (int jj = 0; jj < num_elemental_dofs; jj++)
-        dfe[ii] =
-            dfe[ii] + Ke[ii * num_elemental_dofs + jj] * dirichletValue[ii];
-    for (int ii = 0; ii < num_elemental_dofs; ii++)
-      fe[ii] = fe[ii] + dfe[ii];
-    las->AddToVector(num_elemental_dofs, dof_numbers, &fe[0]);
+    las->AddToVector(num_elemental_dofs, dof_numbers.data(), &fe[0]);
   }
 
-
-  void FEAStep::AssembleIntegratorIntoLAS(LAS * las,
-                                 apf::Field * coordinates)
+  void FEAStep::AssembleIntegratorIntoLAS(LAS * las, apf::Field * coordinates)
 
   {
     if (coordinates == nullptr)
     {
       coordinates = apf_mesh->getCoordinateField();
     }
-    apf::MeshIterator * it = apf_mesh->begin(analysis_dim);
-    apf::MeshEntity * me = nullptr;
-    while ((me = apf_mesh->iterate(it)))
+    apf::MeshIterator * mesh_region_iter = apf_mesh->begin(analysis_dim);
+    while (apf::MeshEntity * mesh_entity = apf_mesh->iterate(mesh_region_iter))
     {
-      if (!apf_mesh->isOwned(me))
+      if (!apf_mesh->isOwned(mesh_entity))
       {
         continue;
       }
-      apf::MeshElement * mlm = apf::createMeshElement(coordinates, me);
-      auto* sys = getIntegrator(me, 0);
-      apf::Element * elm = apf::createElement(sys->getField(), mlm);
+      apf::MeshElement * mlm = apf::createMeshElement(coordinates, mesh_entity);
+      auto * sys = getIntegrator(mesh_entity, 0);
       sys->process(mlm);
-      apf::NewArray<apf::Vector3> dofs;
-      apf::getVectorNodes(elm, dofs);
-      apf::NewArray<int> ids;
-      apf::getElementNumbers(apf_primary_numbering, me, ids);
-      AssembleDOFs(las, sys->numElementalDOFs(), &ids[0], &dofs[0],
-                   &sys->getKe()(0, 0), &sys->getfe()(0),
-                   sys->includesBodyForces(), analysis_dim);
-      apf::destroyElement(elm);
+
+      AssembleDOFs(las, sys->getFieldNumbers(), sys->getFieldValues(),
+                   sys->getKe(), sys->getfe(), sys->includesBodyForces(),
+                   apf::countComponents(sys->getField()));
       apf::destroyMeshElement(mlm);
     }
-    apf_mesh->end(it);
+    apf_mesh->end(mesh_region_iter);
   }
 
 }  // namespace amsi

src/mumfim/analysis/amsiLAS.h

@@ -15,13 +15,17 @@ namespace amsi
     /// add an element val at given row and column in the system
     virtual void AddToMatrix(int row, int col, double value ) = 0;
     /// add a matrix of values to specific locations in the global system, useful for adding elemental system's contributions
-    virtual void AddToMatrix(int num_row, int * rows,
-                             int num_col, int * cols,
-                             double * values) = 0;
+    virtual void AddToMatrix(int num_row,
+                             const int * rows,
+                             int num_col,
+                             const int * cols,
+                             const double * values) = 0;
     /// add an element in the right hand side
     virtual void AddToVector(int row, double value) = 0;
     /// add multiple elements to the right hand side
-    virtual void AddToVector(int num_rows, int * rows, double * values) = 0;
+    virtual void AddToVector(int num_rows,
+                             const int * rows,
+                             const double * values) = 0;
     /// solve the system
     virtual void solve () = 0;
     /// zero the system

src/mumfim/analysis/amsiLinearElasticConstitutive.cc

@@ -30,8 +30,8 @@ namespace amsi
     result(5,5) = mu;
     return result;
   }
-  LinearElasticIntegrator::LinearElasticIntegrator(apf::Field * field, int o, double E, double v)
-    : ElementalSystem(field,o)
+  LinearElasticIntegrator::LinearElasticIntegrator(apf::Field * field, apf::Numbering* numbering, int o, double E, double v)
+    : ElementalSystem(field, numbering, o)
     , C(isotropicLinearElasticityTensor(E,v))
   {
     fs = apf::getShape(f);

src/mumfim/analysis/amsiLinearElasticConstitutive.h

@@ -8,6 +8,7 @@ namespace amsi
   {
   public:
     LinearElasticIntegrator(apf::Field * field,
+                            apf::Numbering* numbering,
                             int o,
                             double E,
                             double v);

src/mumfim/analysis/amsiPETScLAS.cc

@@ -79,7 +79,10 @@ namespace amsi
   {
     MatSetValues(A,1,&row,1,&col,&value,ADD_VALUES);
   }
-  void PetscLAS::AddToMatrix(int num_rows, int * rows, int num_cols, int * cols, double * values)
+  void PetscLAS::AddToMatrix(int num_rows,
+                             const int * rows, int num_cols,
+                             const int * cols,
+                             const double * values)
   {
     MatSetValues(A,num_rows,rows,num_cols,cols,values,ADD_VALUES);
   }
@@ -102,15 +105,17 @@ namespace amsi
     VecSetValue(b_i,row,static_cast<PetscScalar>(value),ADD_VALUES);
     b_assembled = false;
   }
-  void PetscLAS::AddToVector(int num_rows, int * rows, double * values)
+  void PetscLAS::AddToVector(int num_rows,
+                             const int * rows,
+                             const double * values)
   {
     if(!b_addMode)
     {
       VecAssemblyBegin(b_i);
       VecAssemblyEnd(b_i);
       b_addMode = true;
     }
-    VecSetValues(b_i,num_rows,rows,static_cast<PetscScalar*>(values),ADD_VALUES);
+    VecSetValues(b_i,num_rows,rows,static_cast<const PetscScalar*>(values),ADD_VALUES);
     b_assembled = false;
   }
   /**

src/mumfim/analysis/amsiPETScLAS.h

@@ -16,9 +16,9 @@ namespace amsi
     void Reinitialize(int, int, int, int*);
     void Reinitialize(int, int, int);
     void AddToMatrix(int,int,double);
-    void AddToMatrix(int,int*,int,int*,double*);
+    void AddToMatrix(int, const int *,int, const int *, const double *);
     void AddToVector(int,double);
-    void AddToVector(int,int*,double*);
+    void AddToVector(int, const int *, const double *);
     void solve();
     bool Zero();
     bool ZeroMatrix();