17 {
18
20
21 try {
22
23 moab::Core mb_instance;
24 moab::Interface &moab = mb_instance;
25 int rank;
26 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
27
28 PetscBool flg = PETSC_TRUE;
32 if (flg != PETSC_TRUE) {
33 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
34 }
35
36 const char *option;
37 option = "";
39
42
43
45 bit_level0.set(0);
47 CHKERR moab.create_meshset(MESHSET_SET, meshset_level0);
49 0, 3, bit_level0);
52
53
55 3);
56
57
59
60
62 "SPATIAL_POSITION");
64 "SPATIAL_POSITION");
66 "SPATIAL_POSITION");
67
68
70
71
73 "ELASTIC");
74
75
77 bit_level0);
78
79
81
82
85
86
89 &flg);
90 if (flg != PETSC_TRUE) {
92 }
97
100 "SPATIAL_POSITION");
102 "SPATIAL_POSITION");
104 "SPATIAL_POSITION");
106 "NEUAMNN_FE");
108 it)) {
110 CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris,
true);
112 "NEUAMNN_FE");
113 }
117 CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris,
true);
119 "NEUAMNN_FE");
120 }
121
122
124 3);
126 int order_velocity = 1;
128 order_velocity);
130 order_velocity);
132 order_velocity);
134
146 order_velocity);
148 order_velocity);
150 order_velocity);
152
154 CHKERR inertia.setBlocks();
155 CHKERR inertia.addConvectiveMassElement(
"MASS_ELEMENT",
"SPATIAL_VELOCITY",
156 "SPATIAL_POSITION");
157 CHKERR inertia.addVelocityElement(
"VELOCITY_ELEMENT",
"SPATIAL_VELOCITY",
158 "SPATIAL_POSITION");
160 "MASS_ELEMENT");
162 "VELOCITY_ELEMENT");
163
164
166
167 double scale_positions = 2;
168 {
170 double coords[3];
172 dof_ptr)) {
173 if (dof_ptr->get()->getEntType() != MBVERTEX)
174 continue;
176 int dof_rank = dof_ptr->get()->getDofCoeffIdx();
177 double &fval = dof_ptr->get()->getFieldData();
178 if (node != ent) {
179 CHKERR moab.get_coords(&ent, 1, coords);
180 node = ent;
181 }
182 fval = scale_positions * coords[dof_rank];
183 }
184 }
185
186 double scale_velocities = 4;
187 {
189 double coords[3];
191 dof_ptr)) {
192 if (dof_ptr->get()->getEntType() != MBVERTEX)
193 continue;
195 int dof_rank = dof_ptr->get()->getDofCoeffIdx();
196 double &fval = dof_ptr->get()->getFieldData();
197 if (node != ent) {
198 CHKERR moab.get_coords(&ent, 1, coords);
199 node = ent;
200 }
201 fval = scale_velocities * coords[dof_rank];
202 }
203 }
204
205
207
208
210
211
212
215
217
218
222
223
226 "ELASTIC_MECHANICS",
COL, &F);
229 Mat Aij;
231 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("ELASTIC_MECHANICS",
232 &Aij);
233
234 CHKERR inertia.setConvectiveMassOperators(
"SPATIAL_VELOCITY",
235 "SPATIAL_POSITION");
236 CHKERR inertia.setVelocityOperators(
"SPATIAL_VELOCITY",
"SPATIAL_POSITION");
237
238 inertia.getLoopFeMassRhs().ts_F = F;
239 inertia.getLoopFeMassRhs().ts_a = 1;
240 inertia.getLoopFeMassLhs().ts_B = Aij;
241 inertia.getLoopFeMassLhs().ts_a = 1;
242
243 inertia.getLoopFeVelRhs().ts_F = F;
244 inertia.getLoopFeVelRhs().ts_a = 1;
245 inertia.getLoopFeVelLhs().ts_B = Aij;
246 inertia.getLoopFeVelLhs().ts_a = 1;
247
249 inertia.getLoopFeMassRhs());
251 inertia.getLoopFeVelRhs());
252 CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
253 CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
254 CHKERR VecAssemblyBegin(F);
256
258 inertia.getLoopFeMassLhs());
260 inertia.getLoopFeVelLhs());
261 CHKERR MatAssemblyBegin(Aij, MAT_FINAL_ASSEMBLY);
262 CHKERR MatAssemblyEnd(Aij, MAT_FINAL_ASSEMBLY);
263
264 double sum = 0;
266 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"sum = %9.8e\n", sum);
267 double fnorm;
268 CHKERR VecNorm(F, NORM_2, &fnorm);
269 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"fnorm = %9.8e\n", fnorm);
270
271 double mnorm;
272 CHKERR MatNorm(Aij, NORM_1, &mnorm);
273 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"mnorm = %9.8e\n", mnorm);
274
275 if (fabs(sum - 6.27285463e+00) > 1e-8) {
277 }
278 if (fabs(fnorm - 1.28223353e+00) > 1e-6) {
280 }
281 if (fabs(mnorm - 1.31250000e+00) > 1e-6) {
283 }
284
288 }
290
292
293 return 0;
294}
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ MOFEM_ATOM_TEST_INVALID
#define CHKERR
Inline error check.
#define _IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_(MFIELD, NAME, IT)
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 add_ents_to_finite_element_by_bit_ref(const BitRefLevel &bit, const BitRefLevel &mask, const std::string &name, EntityType type, int verb=DEFAULT_VERBOSITY)=0
add TET entities from given refinement level to finite element database given by name
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 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.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
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
const FTensor::Tensor2< T, Dim, Dim > Vec
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
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)
structure grouping operators and data used for calculation of mass (convective) element \ nonlinear_e...
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.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Matrix manager is used to build and partition problems.
Problem manager is used to build and partition problems.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.