v0.14.0
Loading...
Searching...
No Matches
Functions | Variables
damper_jacobian_test.cpp File Reference

Atom test testing calculation of element residual vectors and tangent matrices. More...

#include <BasicFiniteElements.hpp>

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 

Detailed Description

Atom test testing calculation of element residual vectors and tangent matrices.

Definition in file damper_jacobian_test.cpp.

Function Documentation

◆ main()

int main ( int argc,
char * argv[] )

Definition at line 17 of file damper_jacobian_test.cpp.

17 {
18
19 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
20
21 try {
22
23 moab::Core mb_instance;
24 moab::Interface &moab = mb_instance;
25
26 {
27 PetscBool flg = PETSC_TRUE;
28 char mesh_file_name[255];
29 CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
30 mesh_file_name, 255, &flg);
31 ;
32 if (flg != PETSC_TRUE) {
33 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
34 }
35 const char *option;
36 option = ""; //"PARALLEL=BCAST;";//;DEBUG_IO";
37 CHKERR moab.load_file(mesh_file_name, 0, option);
38 }
39
40 MoFEM::Core core(moab);
41 MoFEM::Interface &m_field = core;
42 BitRefLevel bit_level0;
43 bit_level0.set(0);
44 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
45 0, 3, bit_level0);
46 ;
47
48 // Define fields and finite elements
49 {
50
51 // Set approximation fields
52 {
53 // Seed all mesh entities to MoFEM database, those entities can be
54 // potentially used as finite elements or as entities which carry some
55 // approximation field.
56
57 bool check_if_spatial_field_exist =
58 m_field.check_field("SPATIAL_POSITION");
59 CHKERR m_field.add_field("SPATIAL_POSITION", H1,
60 AINSWORTH_LEGENDRE_BASE, 3, MB_TAG_SPARSE,
61 MF_ZERO);
62
63 // meshset consisting all entities in mesh
64 EntityHandle root_set = moab.get_root_set();
65 // add entities to field (root_mesh, i.e. on all mesh entities fields
66 // are approx.)
67 CHKERR m_field.add_ents_to_field_by_type(root_set, MBTET,
68 "SPATIAL_POSITION");
69
70 PetscBool flg;
71 PetscInt order;
72 CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_order", &order,
73 &flg);
74 ;
75 if (flg != PETSC_TRUE) {
76 order = 2;
77 }
78 if (order < 2) {
79 // SETERRQ()
80 }
81
82 CHKERR m_field.set_field_order(root_set, MBTET, "SPATIAL_POSITION",
83 order);
84 CHKERR m_field.set_field_order(root_set, MBTRI, "SPATIAL_POSITION",
85 order);
86 CHKERR m_field.set_field_order(root_set, MBEDGE, "SPATIAL_POSITION",
87 order);
88 CHKERR m_field.set_field_order(root_set, MBVERTEX, "SPATIAL_POSITION",
89 1);
90
91 CHKERR m_field.build_fields();
92
93 // Sett geometry approximation and initial spatial positions
94 // 10 node tets
95 if (!check_if_spatial_field_exist) {
96 Projection10NodeCoordsOnField ent_method_spatial(m_field,
97 "SPATIAL_POSITION");
98 CHKERR m_field.loop_dofs("SPATIAL_POSITION", ent_method_spatial);
99 }
100 }
101
102 // Set finite elements
103 {
104 CHKERR m_field.add_finite_element("DAMPER_FE", MF_ZERO);
106 "SPATIAL_POSITION");
108 "SPATIAL_POSITION");
109
111 "SPATIAL_POSITION");
112
113 EntityHandle root_set = moab.get_root_set();
114 CHKERR m_field.add_ents_to_finite_element_by_type(root_set, MBTET,
115 "DAMPER_FE");
116
117 // build finite elemnts
119 // build adjacencies
120 CHKERR m_field.build_adjacencies(bit_level0);
121 }
122 }
123
124 // Create damper instance
125 KelvinVoigtDamper damper(m_field);
126 {
127
128 int id = 0;
130 damper.blockMaterialDataMap[id];
131 damper.constitutiveEquationMap.insert(
132 id,
134
135 // Set material parameters
136 CHKERR moab.get_entities_by_type(0, MBTET, material_data.tEts);
137 material_data.gBeta = 1;
138 material_data.vBeta = 0.3;
139
140 KelvinVoigtDamper::CommonData &common_data = damper.commonData;
141 common_data.spatialPositionName = "SPATIAL_POSITION";
142 common_data.spatialPositionNameDot = "DOT_SPATIAL_POSITION";
143
144 for (auto &&fe_ptr : {&damper.feRhs, &damper.feLhs}) {
145 // fe_ptr->getOpPtrVector().push_back(
146 // new OpCalculateVectorFieldValuesDot<3>(
147 // common_data.spatialPositionName,
148 // common_data.dataAtGaussTmpPtr));
149 fe_ptr->getOpPtrVector().push_back(
151 common_data.spatialPositionName,
152 common_data.gradDataAtGaussTmpPtr));
153 fe_ptr->getOpPtrVector().push_back(
155 common_data.spatialPositionName, common_data, false, true, false));
156 fe_ptr->getOpPtrVector().push_back(
158 common_data.spatialPositionName,
159 common_data.gradDataAtGaussTmpPtr));
160 fe_ptr->getOpPtrVector().push_back(
162 common_data.spatialPositionName, common_data, false, true, true));
163 }
164
165 // attach tags for each recorder
166 std::vector<int> tags;
167 tags.push_back(1);
168
170 damper.constitutiveEquationMap.at(id);
171
172 // Right hand side operators
173 damper.feRhs.getOpPtrVector().push_back(new KelvinVoigtDamper::OpJacobian(
174 "SPATIAL_POSITION", tags, ce, damper.commonData, true, false));
175 damper.feRhs.getOpPtrVector().push_back(
176 new KelvinVoigtDamper::OpRhsStress(damper.commonData));
177
178 // Left hand side operators
179 damper.feLhs.getOpPtrVector().push_back(new KelvinVoigtDamper::OpJacobian(
180 "SPATIAL_POSITION", tags, ce, damper.commonData, false, true));
181 damper.feLhs.getOpPtrVector().push_back(
182 new KelvinVoigtDamper::OpLhsdxdx(damper.commonData));
183 }
184
185 // Create dm instance
186 DM dm;
187 DMType dm_name = "DMDAMPER";
188 {
189 CHKERR DMRegister_MoFEM(dm_name);
190 CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
191 CHKERR DMSetType(dm, dm_name);
192 CHKERR DMMoFEMCreateMoFEM(dm, &m_field, dm_name, bit_level0);
193 CHKERR DMSetFromOptions(dm);
194 // add elements to dm
195 CHKERR DMMoFEMAddElement(dm, "DAMPER_FE");
196 CHKERR DMSetUp(dm);
197 }
198
199 // Make calculations
200 Mat M;
201 Vec F, U_t;
202 {
206 CHKERR VecZeroEntries(F);
207 CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
208 CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
209 CHKERR VecZeroEntries(U_t);
210 CHKERR VecGhostUpdateBegin(U_t, INSERT_VALUES, SCATTER_FORWARD);
211 CHKERR VecGhostUpdateEnd(U_t, INSERT_VALUES, SCATTER_FORWARD);
212 CHKERR MatZeroEntries(M);
213 damper.feRhs.ts_F = F; // Set right hand side vector manually
214 damper.feRhs.ts_u_t = U_t;
215 damper.feRhs.ts_ctx = TSMethod::CTX_TSSETIFUNCTION;
216 CHKERR DMoFEMLoopFiniteElements(dm, "DAMPER_FE", &damper.feRhs);
217 damper.feLhs.ts_B = M; // Set matrix M
218 damper.feLhs.ts_a = 1.0; // Set time step parameter
219 damper.feLhs.ts_u_t = U_t;
220 CHKERR DMoFEMLoopFiniteElements(dm, "DAMPER_FE", &damper.feLhs);
221 CHKERR VecAssemblyBegin(F);
222 CHKERR VecAssemblyEnd(F);
223 CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
224 CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
225 }
226
227 // See results
228 {
229
230 PetscViewer viewer;
231 CHKERR PetscViewerASCIIOpen(PETSC_COMM_WORLD, "damper_jacobian_test.txt",
232 &viewer);
233 CHKERR PetscViewerDestroy(&viewer);
234 // MatView(M,PETSC_VIEWER_DRAW_WORLD);
235 // std::string wait;
236 // std::cin >> wait;
237 }
238
239 // Clean and destroy
240 {
241 CHKERR VecDestroy(&F);
242 CHKERR VecDestroy(&U_t);
243 CHKERR MatDestroy(&M);
244 CHKERR DMDestroy(&dm);
245 }
246 }
248
250
251 return 0;
252}
static Index< 'M', 3 > M
static char help[]
#define CATCH_ERRORS
Catch errors.
@ MF_ZERO
Definition definitions.h:98
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
@ H1
continuous field
Definition definitions.h:85
#define CHKERR
Inline error check.
constexpr int order
@ F
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:483
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition DMMoFEM.cpp:118
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition DMMoFEM.cpp:1183
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:572
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition DMMoFEM.cpp:1153
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
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 bool check_field(const std::string &name) const =0
check if field is in database
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
char mesh_file_name[255]
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
double vBeta
Poisson ration spring alpha.
double gBeta
Sheer modulus spring alpha.
Common data for nonlinear_elastic_elem model.
boost::shared_ptr< MatrixDouble > gradDataAtGaussTmpPtr
Assemble internal force vector.
Implementation of Kelvin Voigt Damper.
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.
Get field gradients time derivative at integration pts for scalar filed rank 0, i....
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Projection of edge entities with one mid-node on hierarchical basis.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

Variable Documentation

◆ help

char help[] = "...\n\n"
static

Definition at line 15 of file damper_jacobian_test.cpp.