Refactor fem steps

src/mumfim/analysis/ElementalSystem.cc

@@ -12,15 +12,20 @@ namespace amsi
                                             std::vector<int>& dof_numbers, std::vector<double> & values) {
 
     const auto ncomponents = apf::countComponents(field);
-    values.resize(ncomponents*nenodes);
+
     dof_numbers.resize(ncomponents*nenodes);
+    apf::NewArray<int> new_array_numbering(ncomponents*nenodes);
+    apf::getElementNumbers(numbering, mesh_entity, new_array_numbering);
+
+    values.resize(ncomponents*nenodes);
     apf::NewArray<double> new_array_data(ncomponents*nenodes);
     e->getElementNodeData(new_array_data);
+
     for(int node=0; node<nenodes; ++node) {
       for(int component = 0; component < ncomponents; component++){
         int index = node*ncomponents + component;
+        dof_numbers[index] = new_array_numbering[index];
         values[index] = new_array_data[index];
-        dof_numbers[index] = apf::getNumber(numbering,mesh_entity,node, component);
       }
     }
   }

src/mumfim/analysis/ElementalSystem.h

@@ -28,6 +28,8 @@ namespace amsi
     virtual int numElementalDOFs() {return nedofs;}
     virtual apf::DynamicMatrix& getKe() {return Ke;}
     virtual apf::DynamicVector& getfe() {return fe;}
+    virtual void zeroKe(){Ke.zero();}
+    virtual void zerofe(){fe.zero();}
 
     // TODO have getFieldValues and getFieldNumbers return mdspan
     [[nodiscard]] const std::vector<double>& getFieldValues() { return field_values_; }

src/mumfim/analysis/amsiFEA.cc

@@ -224,6 +224,8 @@ namespace amsi
   void FEAStep::AssembleIntegratorIntoLAS(LAS * las, apf::Field * coordinates)
 
   {
+    las->ZeroVector();
+    las->ZeroMatrix();
     if (coordinates == nullptr)
     {
       coordinates = apf_mesh->getCoordinateField();

src/mumfim/macroscale/FEMAnalysis.cc

@@ -3,6 +3,7 @@
 #include <amsiPETScLAS.h>
 #include <apf.h>
 #include <petscsnes.h>
+#include <petscvec.h>
 
 #include <iostream>
 
@@ -56,8 +57,8 @@ namespace mumfim
       // residuals
       const double * sol;
       VecGetArrayRead(solution, &sol);
+
       an->analysis_step_->UpdateDOFs(sol);
-      VecRestoreArrayRead(solution, &sol);
       an->analysis_step_->ApplyBC_Dirichlet();
       an->analysis_step_->RenumberDOFs();
       int gbl, lcl, off;
@@ -70,10 +71,13 @@ namespace mumfim
       VecAssemblyEnd(petsc_las->GetVector());
       MatAssemblyBegin(petsc_las->GetMatrix(), MAT_FINAL_ASSEMBLY);
       MatAssemblyEnd(petsc_las->GetMatrix(), MAT_FINAL_ASSEMBLY);
+
+      Vec rhs = petsc_las->GetVector(); 
+
+      MatMultAdd(petsc_las->GetMatrix(), solution, rhs, residual);
+
       an->analysis_step_->iter();
-      // instead of doing this we can probably pas las->GetVector() to SNESSetFunction
-      VecCopy(petsc_las->GetVector(), residual);
-      // an->analysis_step_->AcceptDOFs();
+      VecRestoreArrayRead(solution, &sol);
     }
     catch (mumfim_error & e)
     {

src/mumfim/macroscale/LinearHeatConductionStep.cc

@@ -29,12 +29,9 @@ namespace mumfim
     apf_primary_numbering = apf::createNumbering(apf_primary_field);
 
     kappa = apf::createIPField(apf_mesh, "kappa", apf::MATRIX, 1);
-
-    // Future expansion
-    //flux = apf::createIPField(apf_mesh, "fluxes", apf::VECTOR, 1);
-    //sources = apf::createIPField(apf_mesh, "sources", apf::SCALAR, 1);
-    //apf::zeroField(sources);
-    //apf::zeroField(flux);
+    apf::zeroField(kappa);
+    // used to place kappa into IPfield by region
+    std::map<int, apf::Matrix3x3> mappa;
 
     amsi::applyUniqueRegionTags(apf_mesh);
 
@@ -51,38 +48,42 @@ namespace mumfim
           mt::GetCategoryByType(model_traits, "material model");
       const auto * continuum_model =
           mt::GetPrimaryCategoryByType(material_model, "continuum model");
-      const auto * tmp =
+      const auto * kappa_mt =
           //mt::GetCategoryModelTraitByType<mt::MatrixMT>(continuum_model,
           mt::GetCategoryModelTraitByType<mt::ScalarMT>(continuum_model,
                                                         "kappa");
-      auto k = (*tmp)();
+      if (kappa_mt == nullptr)
+      {
+        std::cerr << " \"kappa\" (thermal conductivity) is required for "
+                     "the continuum model.\n";
+        MPI_Abort(AMSI_COMM_WORLD, 1);
+      }
+
+      auto k = (*kappa_mt)();
       // This needs to be unique_ptr-ized, as it currently isn't deleted anywhere.
       apf::Matrix3x3 * kappa_r = new apf::Matrix3x3(
         k  , 0.0, 0.0,
         0.0, k  , 0.0,
         0.0, 0.0, k  );
-
-      /*
-      apf::Matrix3x3 * kappa_r = new apf::Matrix3x3(
-        (*tmp)(0,0), (*tmp)(1,0), (*tmp)(2,0),
-        (*tmp)(0,1), (*tmp)(1,1), (*tmp)(2,1),
-        (*tmp)(0,2), (*tmp)(1,2), (*tmp)(2,2)
-      );
-      */
+      // Save for the kappa IPField assignment
+      mappa[tag] = *kappa_r;
 
       std::cout << "continuum model type: " << continuum_model->GetType()
                 << "\n";
-      if (kappa == nullptr)
-      {
-        std::cerr << " \"kappa\" (thermal conductivity) is required for "
-                     "the continuum model.\n";
-        MPI_Abort(AMSI_COMM_WORLD, 1);
-      }
       constitutives[reinterpret_cast<apf::ModelEntity *>(gent)] =
           std::make_unique<LinearHeatIntegrator>(
               apf_primary_field, apf_primary_numbering, kappa_r);
     }
     gmi_end(gmodel, it);
+
+    apf::MeshEntity *ent;
+    auto *mesh_it = mesh->begin(3);
+    while((ent = mesh->iterate(mesh_it)))
+    {
+      int tag = mesh->getModelTag(mesh->toModel(ent));
+      apf::setMatrix(kappa, ent, 0, mappa.at(tag));
+    }
+
     neumann_bcs.push_back(
         amsi::NeumannBCEntry{ .categories = {"heat_flux"},
                               .mt_name = "x",
@@ -132,91 +133,4 @@ namespace mumfim
   {
     return constitutives[apf_mesh->toModel(mesh_entity)].get();
   }
-  /*
-  void LinearHeatConductionStep::step()
-  {
-    iteration = 0;
-  }
-  */
-  /*
-  void LinearHeatConductionStep::iter()
-  {
-    iteration++;
-  }
-  */
-
-  //void LinearHeatConductionStep::UpdateDOFs(const double * sol) { ; }
-
-  /**
-   * get the Neumann boundary condition value on the specified entity
-   * @param ent model entity to query
-   * @param frc will be filled with the value of the NeumannBCs
-   * This really needs to be reimplemented to be heat flux specific, but I
-   * wasn't willing to just throw away the functionality without replacing it.
-   */
-  /*
-  void LinearHeatConductionStep::getFluxOn(apf::ModelEntity * ent, double * frc)
-  {
-    auto model_dimension = apf_mesh->getModelType(ent);
-    auto model_tag = apf_mesh->getModelTag(ent);
-    auto * associated_traits =
-        problem_definition.associated->Find({model_dimension, model_tag});
-    if (associated_traits == nullptr)
-    {
-      std::cerr << "Attempting to get Neumann BC on [" << model_dimension << ","
-                << model_tag << "].";
-      std::cerr << " No boundary conditions associated with that geometry.\n";
-      exit(1);
-    }
-    for (const auto & neumann_bc : neumann_bcs)
-    {
-      const mt::AssociatedCategoryNode * category_node = associated_traits;
-      for (const auto & category : neumann_bc.categories)
-      {
-        category_node = category_node->FindCategoryByType(category);
-        if (category_node == nullptr)
-        {
-          break;
-        }
-      }
-      if (category_node == nullptr)
-      {
-        std::cerr << "Invalid Neumann BC provided\n";
-        exit(1);
-      }
-      const auto * flux_trait =
-          mt::GetCategoryModelTraitByType(category_node, neumann_bc.mt_name);
-      const auto * const_flux_trait = mt::MTCast<mt::VectorMT>(flux_trait);
-      const auto * time_function_flux =
-          mt::MTCast<mt::VectorFunctionMT<1>>(flux_trait);
-      if (const_flux_trait != nullptr)
-      {
-        for (int i = 0; i < 3; ++i)
-        {
-          frc[i] = (*const_flux_trait)(i);
-        }
-      }
-      else if (time_function_flux != nullptr)
-      {
-        for (int i = 0; i < 3; ++i)
-        {
-          frc[i] = (*time_function_flux)(i, T);
-        }
-      }
-      else
-      {
-        std::cerr << "flux must be either a scalar value, or function of time "
-                     "only.\n";
-        exit(1);
-      }
-    }
-  }
-*/
-/*
-  void LinearHeatConductionStep::recoverSecondaryVariables(int)
-  {
-    ;
-  }
-  */
-
-}  // namespace mumfim
+}

src/mumfim/macroscale/LinearHeatIntegrator.cc

@@ -38,26 +38,6 @@ namespace mumfim
     ;
   }
 
-  LinearHeatIntegrator::LinearHeatIntegrator(apf::Field * temperature,apf::Numbering* numbering, apf::Field * kappa)
-    : amsi::ElementalSystem(temperature,numbering, 1)
-    , T_(temperature)
-    , K_(kappa)
-  {
-    ;
-  }
-
-  /*
-  void LinearHeatIntegrator::inElement(apf::MeshElement *me)
-  {
-    e = apf::createElement(T_, me);
-    nenodes = apf::countNodes(e);
-    nedofs = nenodes;  // Assumes scalar values at nodes
-    Ke.setSize(nenodes, nenodes);
-    Ke.zero();
-    fe.setSize(nenodes);
-    fe.zero();
-  }
-  */
 
   void LinearHeatIntegrator::atPoint(apf::Vector3 const &p, double w, double dV)
   {

src/mumfim/macroscale/LinearHeatIntegrator.h

@@ -7,7 +7,6 @@ namespace mumfim{
   class LinearHeatIntegrator : public amsi::ElementalSystem
   {
   apf::Field * T_;
-  apf::Field * K_;
   apf::Field * f_;
   apf::DynamicMatrix D;  // elemental kappa
   apf::Matrix3x3 * getIsotropicKappa(double kii);