v0.14.0
Loading...
Searching...
No Matches
nonlinear_elastic_test_jacobian.cpp
Go to the documentation of this file.
2using namespace MoFEM;
3
4#include <Hooke.hpp>
5#include <NeoHookean.hpp>
6
7namespace bio = boost::iostreams;
8using bio::stream;
9using bio::tee_device;
10
11static char help[] = "...\n\n";
12
13struct OpCheck
18 "SPATIAL_POSITION", UserDataOperator::OPROW),
19 commonData(common_data) {}
20
23
27 const int nb_dofs = data.getFieldData().size();
28 // const int nb_base_functions = data.getN().size2();
29 if (nb_dofs == 0) {
31 }
32 const double eps = 1e-12;
33 const int nb_gauss_pts = data.getN().size1();
34 for (int gg = 0; gg != nb_gauss_pts; gg++) {
37 str0(0, 0), str0(0, 1), str0(0, 2), str0(1, 0), str0(1, 1),
38 str0(1, 2), str0(2, 0), str0(2, 1), str0(2, 2));
40 FTensor::Tensor2<double, 3, 3> t_stress_1(str1[0], str1[1], str1[2],
41 str1[3], str1[4], str1[5],
42 str1[6], str1[7], str1[8]);
43 t_stress_0(i, j) -= t_stress_1(i, j);
44 double nrm2 = t_stress_0(i, j) * t_stress_0(i, j);
45 // PetscPrintf(PETSC_COMM_WORLD,"nrm2 = %3.2e\n",nrm2);
46 if (nrm2 > eps) {
47 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
48 "derivative of energy and stress inconsistency nrm2 = %6.4e",
49 nrm2);
50 }
51 }
53 }
54};
55
56int main(int argc, char *argv[]) {
57
58 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
59
60 try {
61
62 enum materials { HOOKE, NEOHOOKEAN, LASTOP };
63
64 const char *materials_list[] = {"HOOKE", "NEOHOOKEAN"};
65
66 PetscBool flg_test_mat;
67 PetscInt choise_value = HOOKE;
68 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-mat", materials_list,
69 LASTOP, &choise_value, &flg_test_mat);
70
71 moab::Core mb_instance;
72 moab::Interface &moab = mb_instance;
73 int rank;
74 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
75
76 PetscBool flg = PETSC_TRUE;
77 char mesh_file_name[255];
78 CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
79 mesh_file_name, 255, &flg);
80
81 const char *option;
82 option = ""; //"PARALLEL=BCAST;";//;DEBUG_IO";
83 CHKERR moab.load_file(mesh_file_name, 0, option);
84
85 MoFEM::Core core(moab);
86 MoFEM::Interface &m_field = core;
87
88 // ref meshset ref level 0
89 BitRefLevel bit_level0;
90 bit_level0.set(0);
91 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
92 0, 3, bit_level0);
93
94 // Fields
95 CHKERR m_field.add_field("SPATIAL_POSITION", H1, AINSWORTH_LEGENDRE_BASE,
96 3);
97 // add entitities (by tets) to the field
98 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "SPATIAL_POSITION");
99 // set app. order
100 PetscInt order;
101 CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_order", &order,
102 &flg);
103 if (flg != PETSC_TRUE) {
104 order = 1;
105 }
106 CHKERR m_field.set_field_order(0, MBTET, "SPATIAL_POSITION", order);
107 CHKERR m_field.set_field_order(0, MBTRI, "SPATIAL_POSITION", order);
108 CHKERR m_field.set_field_order(0, MBEDGE, "SPATIAL_POSITION", order);
109 CHKERR m_field.set_field_order(0, MBVERTEX, "SPATIAL_POSITION", 1);
110
111 // build field
112 CHKERR m_field.build_fields();
113
114 // use this to apply some strain field to the body (testing only)
115 double scale_positions = 2;
116 {
117 EntityHandle node = 0;
118 double coords[3];
119 for (_IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_(m_field, "SPATIAL_POSITION",
120 dof_ptr)) {
121 if (dof_ptr->get()->getEntType() != MBVERTEX)
122 continue;
123 EntityHandle ent = dof_ptr->get()->getEnt();
124 int dof_rank = dof_ptr->get()->getDofCoeffIdx();
125 double &fval = dof_ptr->get()->getFieldData();
126 if (node != ent) {
127 CHKERR moab.get_coords(&ent, 1, coords);
128 node = ent;
129 }
130 fval = scale_positions * coords[dof_rank];
131 }
132 }
133
134 NonlinearElasticElement elastic(m_field, 1);
135
136 CHKERR elastic.setBlocks(
137 boost::make_shared<NonlinearElasticElement::
138 FunctionsToCalculatePiolaKirchhoffI<double>>(),
139 boost::make_shared<NonlinearElasticElement::
140 FunctionsToCalculatePiolaKirchhoffI<adouble>>());
141 elastic.commonData.spatialPositions = "SPATIAL_POSITION";
142 elastic.commonData.meshPositions = "MESH_NODE_POSITIONS";
143
144 CHKERR elastic.addElement("ELASTIC", "SPATIAL_POSITION");
145
146 // build finite elemnts
148 // build adjacencies
149 CHKERR m_field.build_adjacencies(bit_level0);
150
151 // define problems
152 CHKERR m_field.add_problem("ELASTIC_MECHANICS");
153 // set refinement level for problem
154 CHKERR m_field.modify_problem_ref_level_add_bit("ELASTIC_MECHANICS",
155 bit_level0);
156 // set finite elements for problems
157 CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
158 "ELASTIC");
159
160 ProblemsManager *prb_mng_ptr;
161 CHKERR m_field.getInterface(prb_mng_ptr);
162 // build problem
163 CHKERR prb_mng_ptr->buildProblem("ELASTIC_MECHANICS", true);
164 // partition
165 CHKERR prb_mng_ptr->partitionProblem("ELASTIC_MECHANICS");
166 CHKERR prb_mng_ptr->partitionFiniteElements("ELASTIC_MECHANICS");
167 CHKERR prb_mng_ptr->partitionGhostDofs("ELASTIC_MECHANICS");
168
169 // test nonlinear element
170 {
171 elastic.getLoopFeRhs().getOpPtrVector().push_back(
173 "SPATIAL_POSITION", elastic.commonData));
174 const int tag0 = 1;
175 const int tag1 = 2;
176 std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
177 elastic.setOfBlocks.begin();
178 for (; sit != elastic.setOfBlocks.end(); sit++) {
179 elastic.getLoopFeRhs().getOpPtrVector().push_back(
181 "SPATIAL_POSITION", sit->second, elastic.commonData, tag0,
182 false, false, false));
183 elastic.getLoopFeRhs().getOpPtrVector().push_back(
185 "SPATIAL_POSITION", sit->second, elastic.commonData, tag1, true,
186 false, false, false));
187 }
188 // Run opperator to check consistency
189 elastic.getLoopFeRhs().getOpPtrVector().push_back(
190 new OpCheck(elastic.commonData));
191
192 CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "ELASTIC",
193 elastic.getLoopFeRhs());
194 }
195
196 // test materials
197 if (flg_test_mat == PETSC_TRUE) {
198 PetscPrintf(PETSC_COMM_WORLD, "Testing %s\n",
199 materials_list[choise_value]);
200 std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
201 elastic.setOfBlocks.begin();
202 for (; sit != elastic.setOfBlocks.end(); sit++) {
203 switch (choise_value) {
204 case HOOKE:
205 sit->second.materialDoublePtr = boost::make_shared<Hooke<double>>();
206 sit->second.materialAdoublePtr = boost::make_shared<Hooke<adouble>>();
207 break;
208 case NEOHOOKEAN:
209 sit->second.materialDoublePtr =
210 boost::make_shared<NeoHookean<double>>();
211 sit->second.materialAdoublePtr =
212 boost::make_shared<NeoHookean<adouble>>();
213 }
214 }
215 }
216 }
218
220
221 return 0;
222}
Implementation of linear elastic material.
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Implementation of Neo-Hookean elastic material.
int main()
Definition: adol-c_atom.cpp:46
static const double eps
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr int order
#define _IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_(MFIELD, NAME, IT)
Definition: Interface.hpp:1843
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
MoFEMErrorCode partitionProblem(const std::string name, int verb=VERBOSE)
partition problem dofs (collective)
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elements
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
char mesh_file_name[255]
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
static char help[]
Managing BitRefLevels.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
@ OPROW
operator doWork function is executed on FE rows
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Problem manager is used to build and partition problems.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
common data used by volume elements
std::vector< MatrixDouble3by3 > sTress
Calculate explicit derivative of free energy.
structure grouping operators and data used for calculation of nonlinear elastic element
MoFEMErrorCode addElement(const std::string element_name, const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false)
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
MoFEMErrorCode setBlocks(boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< double > > materialDoublePtr, boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< adouble > > materialAdoublePtr)
MyVolumeFE & getLoopFeRhs()
get rhs volume element
OpCheck(NonlinearElasticElement::CommonData &common_data)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
NonlinearElasticElement::CommonData & commonData
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.