diff --git a/src/mumfim/analysis/ElementalSystem.cc b/src/mumfim/analysis/ElementalSystem.cc index d35edd97..724a85e5 100644 --- a/src/mumfim/analysis/ElementalSystem.cc +++ b/src/mumfim/analysis/ElementalSystem.cc @@ -1,38 +1,66 @@ #include "ElementalSystem.h" +#include +#include + 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& dof_numbers, std::vector & 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 #include +#include 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 field_values_; + std::vector 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& getFieldValues() { return field_values_; } + [[ nodiscard ]] const std::vector & getFieldNumbers() { return field_numbers_; } virtual apf::Field * getField() {return f;} }; } diff --git a/src/mumfim/analysis/amsiFEA.cc b/src/mumfim/analysis/amsiFEA.cc index 488619fe..85d964e8 100644 --- a/src/mumfim/analysis/amsiFEA.cc +++ b/src/mumfim/analysis/amsiFEA.cc @@ -1,5 +1,6 @@ #include "amsiFEA.h" +#include #include #include #include @@ -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(*component)}; + return {.associated = mt::AssociateModel(component), + .unassociated = std::make_shared(*component)}; } FEAStep::FEAStep(apf::Mesh * mesh, @@ -187,114 +190,61 @@ namespace amsi offset = first_local_dof; } - template - 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 & dof_numbers, + const std::vector & 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(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 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 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 dofs; - apf::getVectorNodes(elm, dofs); - apf::NewArray 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 diff --git a/src/mumfim/analysis/amsiLAS.h b/src/mumfim/analysis/amsiLAS.h index 3df67655..4da59261 100644 --- a/src/mumfim/analysis/amsiLAS.h +++ b/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 diff --git a/src/mumfim/analysis/amsiLinearElasticConstitutive.cc b/src/mumfim/analysis/amsiLinearElasticConstitutive.cc index ad3a8b52..1fdbcd42 100644 --- a/src/mumfim/analysis/amsiLinearElasticConstitutive.cc +++ b/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); diff --git a/src/mumfim/analysis/amsiLinearElasticConstitutive.h b/src/mumfim/analysis/amsiLinearElasticConstitutive.h index b53b906f..5ecf721f 100644 --- a/src/mumfim/analysis/amsiLinearElasticConstitutive.h +++ b/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); diff --git a/src/mumfim/analysis/amsiPETScLAS.cc b/src/mumfim/analysis/amsiPETScLAS.cc index 2ab6cf10..3e6d6f00 100644 --- a/src/mumfim/analysis/amsiPETScLAS.cc +++ b/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,7 +105,9 @@ namespace amsi VecSetValue(b_i,row,static_cast(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) { @@ -110,7 +115,7 @@ namespace amsi VecAssemblyEnd(b_i); b_addMode = true; } - VecSetValues(b_i,num_rows,rows,static_cast(values),ADD_VALUES); + VecSetValues(b_i,num_rows,rows,static_cast(values),ADD_VALUES); b_assembled = false; } /** diff --git a/src/mumfim/analysis/amsiPETScLAS.h b/src/mumfim/analysis/amsiPETScLAS.h index d3610d27..8fb2ce67 100644 --- a/src/mumfim/analysis/amsiPETScLAS.h +++ b/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(); diff --git a/src/mumfim/analysis/apfMeasure.cc b/src/mumfim/analysis/apfMeasure.cc index 9bf2891d..7db81812 100644 --- a/src/mumfim/analysis/apfMeasure.cc +++ b/src/mumfim/analysis/apfMeasure.cc @@ -59,8 +59,8 @@ namespace amsi { double vol; public: - MeasureDisplaced(apf::Field* field, int o) - : ElementalSystem(field, o) + MeasureDisplaced(apf::Field* field, apf::Numbering* numbering, int o) + : ElementalSystem(field, numbering, o) , dim(apf::getMesh(field)->getDimension()) , ent_dim(-1) , fs(NULL) @@ -119,8 +119,8 @@ namespace amsi { }; class MeasureDisplacedFromSurf : public amsi::ElementalSystem { public: - MeasureDisplacedFromSurf(apf::Field* field, int o, int normal_dir) - : ElementalSystem(field, o), dim(0), vol(0), norm_dir(normal_dir) + MeasureDisplacedFromSurf(apf::Field* field, apf::Numbering* numbering, int o, int normal_dir) + : ElementalSystem(field, numbering, o), dim(0), vol(0), norm_dir(normal_dir) { } void inElement(apf::MeshElement* me) @@ -196,10 +196,13 @@ namespace amsi { vol = amsi::comm_sum(vol); return vol; } - double measureDisplacedModelEntity(apf::ModelEntity* mdl_ent, apf::Field* u) + + double measureDisplacedModelEntity(apf::ModelEntity * mdl_ent, + apf::Field * u, + apf::Numbering * numbering) { double vol = 0.0; - MeasureDisplaced elvol(u, 1); + MeasureDisplaced elvol(u, numbering, 1); apf::Mesh* msh = apf::getMesh(u); int dim = msh->getModelType(mdl_ent); apf::MeshEntity* ent; @@ -217,10 +220,13 @@ namespace amsi { vol = amsi::comm_sum(vol); return vol; } - double measureDisplacedMeshEntity(apf::MeshEntity* ent, apf::Field* u) + + double measureDisplacedMeshEntity(apf::MeshEntity * ent, + apf::Field * u, + apf::Numbering * numbering) { // todo : derive integration order from field order - MeasureDisplaced elemental_volume(u, 1); + MeasureDisplaced elemental_volume(u, numbering, 1); apf::Mesh* msh = apf::getMesh(u); apf::MeshElement* mlm = apf::createMeshElement(msh, ent); elemental_volume.process(mlm); @@ -228,11 +234,14 @@ namespace amsi { apf::destroyMeshElement(mlm); return vol; } - double measureDisplacedMeshEntity_greens(apf::MeshEntity* ent, apf::Field* u, - int norm_dir) + + double measureDisplacedMeshEntity_greens(apf::MeshEntity * ent, + apf::Field * u, + int norm_dir, + apf::Numbering * numbering) { // todo : derive integration order from field order - MeasureDisplacedFromSurf elemental_volume(u, 1, norm_dir); + MeasureDisplacedFromSurf elemental_volume(u, numbering, 1, norm_dir); apf::Mesh* msh = apf::getMesh(u); apf::MeshElement* mlm = apf::createMeshElement(msh, ent); elemental_volume.process(mlm); @@ -254,7 +263,7 @@ namespace amsi { if (to_compare == mdl_ent) { int sd = side(mdl_ent, msh, ent); // assert(sd == 1); - vol += measureDisplacedMeshEntity_greens(ent, u, sd); + vol += measureDisplacedMeshEntity_greens(ent, u, sd, nullptr); } } } diff --git a/src/mumfim/analysis/apfMeasure.h b/src/mumfim/analysis/apfMeasure.h index 0328703b..6cb28fec 100644 --- a/src/mumfim/analysis/apfMeasure.h +++ b/src/mumfim/analysis/apfMeasure.h @@ -10,14 +10,18 @@ namespace amsi void faceNormal(const apf::Vector3 & pt_a, const apf::Vector3 & pt_b, const apf::Vector3 & pt_c, apf::Vector3 & n); int side(apf::ModelEntity * mdl_ent, apf::Mesh * msh, apf::MeshEntity * fc); double measureModelEntity(apf::ModelEntity * mdl_ent, apf::Mesh * msh); - double measureDisplacedModelEntity(apf::ModelEntity * mdl_ent, apf::Field * u); + double measureDisplacedModelEntity(apf::ModelEntity * mdl_ent, + apf::Field * u, + apf::Numbering * numbering); /** * @brief Measure the mesh entity w.r.t the entity's dimension, length, area, or volume * @param ent The MeshEntity to measure * @param u The displacement Field to apply to ent * @note Only implemented for regions which gives volume at the moment */ - double measureDisplacedMeshEntity(apf::MeshEntity * ent, apf::Field * u); + double measureDisplacedMeshEntity(apf::MeshEntity * ent, + apf::Field * u, + apf::Numbering * numbering); double measureDisplacedModelEntity_greens(apf::ModelEntity * mdl_ent, apf::Field * u); template double measureDisplacedModelEntities(I b, I e, apf::Field * u); diff --git a/src/mumfim/analysis/apfMeasure_impl.h b/src/mumfim/analysis/apfMeasure_impl.h index aa268e15..fe8f449a 100644 --- a/src/mumfim/analysis/apfMeasure_impl.h +++ b/src/mumfim/analysis/apfMeasure_impl.h @@ -7,7 +7,7 @@ namespace amsi { double v = 0.0; for(auto i = b; i != e; ++i) - v += measureDisplacedModelEntity(*i,u); + v += measureDisplacedModelEntity(*i, u, nullptr); return v; } } diff --git a/src/mumfim/macroscale/HolmesMowIntegrator.h b/src/mumfim/macroscale/HolmesMowIntegrator.h index 545b5283..eb145b7d 100644 --- a/src/mumfim/macroscale/HolmesMowIntegrator.h +++ b/src/mumfim/macroscale/HolmesMowIntegrator.h @@ -14,10 +14,11 @@ namespace mumfim public: HolmesMowIntegrator(NonlinearTissueStep * n, apf::Field * field, + apf::Numbering* numbering, double shear_modulus, double poisson_ratio, int o) - : ElementalSystem(field,o) + : ElementalSystem(field, numbering, o) , current_integration_point(0) , analysis(n) , dim(0) diff --git a/src/mumfim/macroscale/LinearHeatConductionStep.cc b/src/mumfim/macroscale/LinearHeatConductionStep.cc index 4dd544ae..b99d4606 100644 --- a/src/mumfim/macroscale/LinearHeatConductionStep.cc +++ b/src/mumfim/macroscale/LinearHeatConductionStep.cc @@ -79,7 +79,8 @@ namespace mumfim MPI_Abort(AMSI_COMM_WORLD, 1); } constitutives[reinterpret_cast(gent)] = - std::make_unique(apf_primary_field, kappa_r); + std::make_unique( + apf_primary_field, apf_primary_numbering, kappa_r); } gmi_end(gmodel, it); neumann_bcs.push_back( diff --git a/src/mumfim/macroscale/LinearHeatIntegrator.cc b/src/mumfim/macroscale/LinearHeatIntegrator.cc index c5883a38..145fe3cb 100644 --- a/src/mumfim/macroscale/LinearHeatIntegrator.cc +++ b/src/mumfim/macroscale/LinearHeatIntegrator.cc @@ -30,16 +30,16 @@ namespace mumfim return K; } - LinearHeatIntegrator::LinearHeatIntegrator(apf::Field * temperature, apf::Matrix3x3 * kappa) - : amsi::ElementalSystem(temperature,1) + LinearHeatIntegrator::LinearHeatIntegrator(apf::Field * temperature, apf::Numbering* numbering, apf::Matrix3x3 * kappa) + : amsi::ElementalSystem(temperature,numbering, 1) , T_(temperature) , D(apf::fromMatrix(*kappa)) { ; } - LinearHeatIntegrator::LinearHeatIntegrator(apf::Field * temperature, apf::Field * kappa) - : amsi::ElementalSystem(temperature,1) + LinearHeatIntegrator::LinearHeatIntegrator(apf::Field * temperature,apf::Numbering* numbering, apf::Field * kappa) + : amsi::ElementalSystem(temperature,numbering, 1) , T_(temperature) , K_(kappa) { @@ -62,7 +62,6 @@ namespace mumfim // Using TJHughes, pp 68-70, add nonhomogenous part later // fe holds the residual vector, because of the incremental formulation apf::NewArray Na; apf::getShapeValues(e, p, Na); - apf::NewArray nodalTheta; apf::getScalarNodes(e, nodalTheta); apf::NewArray Ba; apf::getShapeGrads(e, p, Ba); apf::DynamicMatrix B(3, nenodes); @@ -74,7 +73,7 @@ namespace mumfim B(0, i) = BT(i, 0) = Ba[i][0]; B(1, i) = BT(i, 1) = Ba[i][1]; B(2, i) = BT(i, 2) = Ba[i][2]; - Theta(i) = nodalTheta[i] * Na[i]; + Theta(i) = field_values_[i] * Na[i]; } apf::DynamicMatrix BT_D_B(nenodes,nenodes); diff --git a/src/mumfim/macroscale/LinearHeatIntegrator.h b/src/mumfim/macroscale/LinearHeatIntegrator.h index e9dc403a..19839782 100644 --- a/src/mumfim/macroscale/LinearHeatIntegrator.h +++ b/src/mumfim/macroscale/LinearHeatIntegrator.h @@ -14,8 +14,8 @@ namespace mumfim{ apf::Matrix3x3 * getFullKappa(double xx, double yy, double zz, double xy, double yz, double zx); public: - LinearHeatIntegrator(apf::Field *T, apf::Field * K); - LinearHeatIntegrator(apf::Field *T, apf::Matrix3x3 * K_r); // Piecewise constant Kappa + LinearHeatIntegrator(apf::Field *T, apf::Numbering* numbering, apf::Field * K); + LinearHeatIntegrator(apf::Field *T, apf::Numbering* numbering, apf::Matrix3x3 * K_r); // Piecewise constant Kappa void inElement(apf::MeshElement *me); bool includesBodyForces() final { return true; } void atPoint(apf::Vector3 const &p, diff --git a/src/mumfim/macroscale/LinearTissueStep.cc b/src/mumfim/macroscale/LinearTissueStep.cc index 9d1c5fbd..6f84eb40 100644 --- a/src/mumfim/macroscale/LinearTissueStep.cc +++ b/src/mumfim/macroscale/LinearTissueStep.cc @@ -85,7 +85,7 @@ namespace mumfim exit(1); } constitutives[tag] = std::make_unique( - apf_primary_field, 1, (*youngs_modulus)(), (*poisson_ratio)()); + apf_primary_field, apf_primary_numbering, 1, (*youngs_modulus)(), (*poisson_ratio)()); } gmi_end(gmodel, it); } diff --git a/src/mumfim/macroscale/MultiscaleTissueStep.cc b/src/mumfim/macroscale/MultiscaleTissueStep.cc index b236f4ac..edb489cb 100644 --- a/src/mumfim/macroscale/MultiscaleTissueStep.cc +++ b/src/mumfim/macroscale/MultiscaleTissueStep.cc @@ -116,7 +116,7 @@ namespace mumfim ornt_3D = apf::createIPField(apf_mesh, "ornt_tens_3D", apf::MATRIX, 1); ornt_2D = apf::createIPField(apf_mesh, "ornt_tens_2D", apf::MATRIX, 1); mltscl = new ULMultiscaleIntegrator(&fo_cplg, microscale_strain, microscale_stress, apf_primary_field, - microscale_dfm_grd, 1,minimum_stiffness); + microscale_dfm_grd, apf_primary_numbering, 1,minimum_stiffness); M2m_id = amsi::getRelationID(multiscale_.getMultiscaleManager(), multiscale_.getScaleManager(), "macro", "micro_fo"); diff --git a/src/mumfim/macroscale/NeoHookeanIntegrator.h b/src/mumfim/macroscale/NeoHookeanIntegrator.h index cf709155..3d599325 100644 --- a/src/mumfim/macroscale/NeoHookeanIntegrator.h +++ b/src/mumfim/macroscale/NeoHookeanIntegrator.h @@ -26,10 +26,11 @@ namespace mumfim apf::Field * current_coords, apf::Field * cauchy_stress_field, apf::Field * green_lagrange_strain_field, + apf::Numbering* numbering, double youngs_modulus, double poisson_ratio, int o) - : ElementalSystem(displacements, o) + : ElementalSystem(displacements, numbering, o) , current_integration_point(0) , dim(0) , dfm_grd_fld(dfm_grd) diff --git a/src/mumfim/macroscale/NonlinearTissueStep.cc b/src/mumfim/macroscale/NonlinearTissueStep.cc index 847f63d9..3627909f 100644 --- a/src/mumfim/macroscale/NonlinearTissueStep.cc +++ b/src/mumfim/macroscale/NonlinearTissueStep.cc @@ -126,7 +126,7 @@ namespace mumfim // solve? constitutives[reinterpret_cast(gent)] = std::make_unique( - apf_primary_field, dfm_grd, current_coords, strs, strn, + apf_primary_field, dfm_grd, current_coords, strs, strn, apf_primary_numbering, (*youngs_modulus)(), (*poisson_ratio)(), 1); } else if (continuum_model->GetType() == "transverse_isotropic") diff --git a/src/mumfim/macroscale/ULMultiscaleIntegrator.h b/src/mumfim/macroscale/ULMultiscaleIntegrator.h index 8baf6ba6..1fff7367 100644 --- a/src/mumfim/macroscale/ULMultiscaleIntegrator.h +++ b/src/mumfim/macroscale/ULMultiscaleIntegrator.h @@ -20,9 +20,10 @@ namespace mumfim apf::Field * strs, apf::Field * u, apf::Field * dfm_grd, + apf::Numbering* numbering, int o, double minimum_stiffness) - : ElementalSystem(u, o) + : ElementalSystem(u, numbering, o) , current_integration_point(0) , coupling(r) , strain_field(strn) diff --git a/src/mumfim/macroscale/UpdatedLagrangianMaterialIntegrator.h b/src/mumfim/macroscale/UpdatedLagrangianMaterialIntegrator.h index 90da1d35..533b19a7 100644 --- a/src/mumfim/macroscale/UpdatedLagrangianMaterialIntegrator.h +++ b/src/mumfim/macroscale/UpdatedLagrangianMaterialIntegrator.h @@ -23,8 +23,9 @@ namespace mumfim apf::Field * strs, apf::Field * u, apf::Field * dfm_grd, + apf::Numbering* numbering, int o) - : ElementalSystem(u, o) + : ElementalSystem(u, numbering, o) , current_integration_point(0) , strain_field(strn) , stress_field(strs) diff --git a/test/macroscale/UpdatedLagrangianMaterialIntegrator.cc b/test/macroscale/UpdatedLagrangianMaterialIntegrator.cc index ab41e70a..7274b932 100644 --- a/test/macroscale/UpdatedLagrangianMaterialIntegrator.cc +++ b/test/macroscale/UpdatedLagrangianMaterialIntegrator.cc @@ -3,6 +3,7 @@ #include #include #include "TestSupport.h" +#include using test::compare_dynamic_matrices; using test::compare_dynamic_vectors; @@ -23,6 +24,8 @@ TEST_CASE("Compare UL Integrators") apf::zeroField(strn); auto * dfm_grd = apf::createIPField(box_mesh, "F", apf::MATRIX, 1); apf::zeroField(dfm_grd); + auto* numbering = apf::createNumbering(displacement); + NaiveOrder(numbering); // set the mesh to have some random displacements auto * it = box_mesh->begin(0); while (auto * entity = box_mesh->iterate(it)) @@ -37,15 +40,15 @@ TEST_CASE("Compare UL Integrators") constexpr double poissons_ratio = 0.33; // create Neohookean integrator mumfim::NeoHookeanIntegrator neohookean(displacement, dfm_grd, current_coords, - strs, strn, youngs_modulus, poissons_ratio, - 1); + strs, strn, numbering, youngs_modulus, + poissons_ratio, 1); // UpdatedLagrangean with NeohookeanMaterial mumfim::UpdatedLagrangianMaterialIntegrator updated_lagrangian( [=](apf::Matrix3x3 F, apf::MeshEntity *, int) { auto shear_modulus = youngs_modulus / (2.0 * (1.0 + poissons_ratio)); return mumfim::NeohookeanMaterial(shear_modulus, poissons_ratio, F); }, - strn, strs, displacement, dfm_grd, 1); + strn, strs, displacement, dfm_grd, numbering, 1); // process mesh it = box_mesh->begin(3); while(apf::MeshEntity * ent = box_mesh->iterate(it)) {