v0.13.2
Loading...
Searching...
No Matches
mass_calculation.cpp
Go to the documentation of this file.
1/* This file is part of MoFEM.
2* \ingroup bone_remodelling
3
4* MoFEM is free software: you can redistribute it and/or modify it under
5* the terms of the GNU Lesser General Public License as published by the
6* Free Software Foundation, either version 3 of the License, or (at your
7* option) any later version.
8*
9* MoFEM is distributed in the hope that it will be useful, but WITHOUT
10* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
12* License for more details.
13*
14* You should have received a copy of the GNU Lesser General Public
15* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
16
17#include <MoFEM.hpp>
18using namespace MoFEM;
20
21#ifdef WITH_METAIO
22#include <metaImage.h>
23#endif
24
25#include <moab/Skinner.hpp>
26#include <moab/AdaptiveKDTree.hpp>
27
28#include <cholesky.hpp>
29#include <DensityMaps.hpp>
30using namespace BoneRemodeling;
31
32#include <PostProcOnRefMesh.hpp>
33
34static char help[] = "\n";
35
36
37
38#ifndef WITH_METAIO
39 #error "You need compile users modules with MetaIO -DWITH_METAIO=1"
40#endif
41
42int main(int argc, char *argv[]) {
43 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
44
45 try {
46
47 moab::Core mb_instance;
48 moab::Interface& moab = mb_instance;
49
50 char mesh_file_name[255];
51 PetscBool flg_file = PETSC_FALSE;
52 char meta_file_name[255];
53 PetscBool meta_flg_file = PETSC_FALSE;
54 PetscInt order = 1;
55 PetscBool is_atom_test = PETSC_FALSE;
56 // char approximation_base[255];
57
58 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Bone remodeling","none"); CHKERRQ(ierr);
59 ierr = PetscOptionsString(
60 "-my_file",
61 "mesh file name","",
62 "mesh.h5m",mesh_file_name,255,&flg_file
63 ); CHKERRQ(ierr);
64
65 ierr = PetscOptionsString(
66 "-meta_image",
67 "meta image name","",
68 "file.mhd",meta_file_name,255,&meta_flg_file
69 ); CHKERRQ(ierr);
70
71 // ierr = PetscOptionsString(
72 // "-base",
73 // "approximation base name","",
74 // "",approximation_base,255,PETSC_NULL
75 // ); CHKERRQ(ierr);
76
77 ierr = PetscOptionsBool(
78 "-my_is_atom_test",
79 "is used with testing, exit with error when diverged","",
80 PETSC_FALSE,&is_atom_test,PETSC_NULL
81 ); CHKERRQ(ierr);
82
83 ierr = PetscOptionsInt(
84 "-my_order",
85 "default approximation order","",
86 1,&order,PETSC_NULL
87 ); CHKERRQ(ierr);
88
89 ierr = PetscOptionsEnd(); CHKERRQ(ierr);
90
91 // #ifdef WITH_METAIO
92
93 MetaImage meta_file(meta_file_name);
94 cout << "Quantity " << meta_file.Quantity() << endl;
95 const int *dim_size = meta_file.DimSize();
96 cout << "Image dimension: " << dim_size[0] << " " << dim_size[1] << " " << dim_size[2] << endl;
97 const double *elemet_size = meta_file.ElementSize();
98 cout << "Image element size: " << elemet_size[0] << " " << elemet_size[1] << " " << elemet_size[2] << endl;
99 const double *offset = meta_file.Offset();
100 cout << "Image offset: " << offset[0] << " " << offset[1] << " " << offset[2] << endl;
101 const double * transform_matrix = meta_file.TransformMatrix();
102 cout << "Image Transform Matrix:\n";
103 for(int rr = 0;rr<3;rr++) {
104 for(int cc = 0;cc<3;cc++) {
105 cout << transform_matrix[rr*3+cc] << " ";
106 }
107 cout << endl;
108 }
109
110 //Read parameters from line command
111 if(flg_file != PETSC_TRUE) {
112 SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
113 }
114 const char *option;
115 option = "";
116 rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
117 ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
118 if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
119
120 MoFEM::Core core(moab);
121 MoFEM::Interface& m_field = core;
122 //meshset consisting all entities in mesh
123 EntityHandle root_set = moab.get_root_set();
124 // Seed all mesh entities to MoFEM database, those entities can be potentially used as finite elements
125 // or as entities which carry some approximation field.
126 BitRefLevel bit_level0;
127 bit_level0.set(0);
128 ierr = m_field.seed_ref_level_3D(root_set,bit_level0); CHKERRQ(ierr);
129 // Add field
130 ierr = m_field.add_field("RHO",H1,AINSWORTH_LEGENDRE_BASE,1); CHKERRQ(ierr);
131 ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);
132
133 // Add meshsets with material and boundary conditions
134 auto *meshsets_manager_ptr = m_field.getInterface<MeshsetsManager>();
135 ierr = meshsets_manager_ptr->setMeshsetFromFile(); CHKERRQ(ierr);
136
137 PetscPrintf(PETSC_COMM_WORLD,"Read meshset. Added meshsets for bc.cfg\n");
138 for(_IT_CUBITMESHSETS_FOR_LOOP_(m_field,it)) {
139 PetscPrintf(
140 PETSC_COMM_WORLD,
141 "%s",static_cast<std::ostringstream&>(std::ostringstream().seekp(0) << *it << endl).str().c_str()
142 );
143 cerr << *it << endl;
144 }
145
146 Range tets_all;
147 rval = m_field.get_moab().get_entities_by_type(0,MBTET,tets_all); CHKERRQ_MOAB(rval);
148
150 string name = it->getName();
151 if (name.compare(0,14,"NO_REMODELLING") == 0) {
152 //nOremodellingBlock = true;
153 EntityHandle meshset = it->getMeshset();
154 Range tets_block;
155 rval = m_field.get_moab().get_entities_by_type(
156 meshset,MBTET,tets_block,true
157 ); CHKERRQ_MOAB(rval);
158 tets_all = subtract(tets_all,tets_block);
159 }
160 }
161
162
163 // if ( strcmp(approximation_base, "lobatto") == 0){
164 // ierr = m_field.add_field("RHO",L2,AINSWORTH_LOBATTO_BASE,1); CHKERRQ(ierr);
165 // } else if ( strcmp(approximation_base, "berstein") == 0){ //not implemented yet
166 // ierr = m_field.add_field("RHO",L2,AINSOWRTH_BERNSTEIN_BEZIER_BASE,1); CHKERRQ(ierr);
167 // }
168 // else {
169 // ierr = m_field.add_field("RHO",L2,AINSWORTH_LEGENDRE_BASE,1); CHKERRQ(ierr);
170 // }
171 // Add entities to field (root_mesh, i.e. on all mesh entities fields are approx.)
172 // On those entities and lower dimension entities approximation is spaned,
173 ierr = m_field.add_ents_to_field_by_type(tets_all,MBTET,"RHO"); CHKERRQ(ierr);
174 ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
175
176
177 {
178 ierr = m_field.add_field("DISTANCE_FROM_SURFACE",H1,AINSWORTH_LEGENDRE_BASE,1); CHKERRQ(ierr);
179 Range vertices;
180 rval = moab.get_entities_by_type(root_set,MBVERTEX,vertices); CHKERRQ_MOAB(rval);
181 ierr = m_field.add_ents_to_field_by_type(vertices,MBTET,"DISTANCE_FROM_SURFACE"); CHKERRQ(ierr);
182 Range tets;
183 rval = moab.get_entities_by_type(root_set,MBTET,tets); CHKERRQ_MOAB(rval);
184 Range edges;
185 rval = moab.get_adjacencies(tets,1,true,edges,moab::Interface::UNION); CHKERRQ_MOAB(rval);
186 ierr = m_field.add_ents_to_field_by_type(edges,MBEDGE,"DISTANCE_FROM_SURFACE"); CHKERRQ(ierr);
187 }
188
189 // Set approximation oder
190 ierr = m_field.set_field_order(root_set,MBVERTEX,"RHO",1); CHKERRQ(ierr);
191 ierr = m_field.set_field_order(root_set,MBEDGE,"RHO",order); CHKERRQ(ierr);
192 ierr = m_field.set_field_order(root_set,MBTRI,"RHO",order); CHKERRQ(ierr);
193 ierr = m_field.set_field_order(root_set,MBTET,"RHO",order); CHKERRQ(ierr);
194 ierr = m_field.set_field_order(root_set,MBEDGE,"DISTANCE_FROM_SURFACE",2); CHKERRQ(ierr);
195 ierr = m_field.set_field_order(root_set,MBVERTEX,"DISTANCE_FROM_SURFACE",1); CHKERRQ(ierr);
196
197 // Apply 2nd order geometry approximation only on skin
198 {
199
200 Skinner skin(&moab);
201 Range faces,tets;
202 rval = moab.get_entities_by_type(0,MBTET,tets); CHKERRQ_MOAB(rval);
203 rval = skin.find_skin(0,tets,false,faces); CHKERRQ_MOAB(rval);
204 Range edges;
205 rval = moab.get_adjacencies(
206 faces,1,false,edges,moab::Interface::UNION
207 ); CHKERRQ_MOAB(rval);
208 ierr = m_field.set_field_order(edges,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
209 }
210 ierr = m_field.set_field_order(root_set,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
211
212 ierr = m_field.build_fields(); CHKERRQ(ierr);
213
214 {
215 Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
216 ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
217 }
218
219 SurfaceKDTree surface_distance(m_field);
220 {
221 Range tets;
222 rval = moab.get_entities_by_type(root_set,MBTET,tets); CHKERRQ_MOAB(rval);
223 ierr = surface_distance.takeASkin(tets); CHKERRQ(ierr);
224 ierr = surface_distance.buildTree(); CHKERRQ(ierr);
225 ierr = surface_distance.setDistanceFromSurface("DISTANCE_FROM_SURFACE"); CHKERRQ(ierr);
226 }
227
228 // Add finite element (this defines element declaration comes later)
229 ierr = m_field.add_finite_element("DENSITY_MAP_ELEMENT"); CHKERRQ(ierr);
230 ierr = m_field.modify_finite_element_add_field_row("DENSITY_MAP_ELEMENT","RHO"); CHKERRQ(ierr);
231 ierr = m_field.modify_finite_element_add_field_col("DENSITY_MAP_ELEMENT","RHO"); CHKERRQ(ierr);
232 ierr = m_field.modify_finite_element_add_field_data("DENSITY_MAP_ELEMENT","RHO"); CHKERRQ(ierr);
233 ierr = m_field.modify_finite_element_add_field_data("DENSITY_MAP_ELEMENT","DISTANCE_FROM_SURFACE"); CHKERRQ(ierr);
234 ierr = m_field.modify_finite_element_add_field_data("DENSITY_MAP_ELEMENT","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
235
236 // Add entities to that element
237 ierr = m_field.add_ents_to_finite_element_by_type(tets_all,MBTET,"DENSITY_MAP_ELEMENT"); CHKERRQ(ierr);
238
239 //build finite elements
240 ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
241 //build adjacencies between elements and degrees of freedom
242 ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
243
244 DM dm_rho;
245 DMType dm_name = "DM_RHO";
246 // Register DM problem
247 ierr = DMRegister_MoFEM(dm_name); CHKERRQ(ierr);
248 ierr = DMCreate(PETSC_COMM_WORLD,&dm_rho);CHKERRQ(ierr);
249 ierr = DMSetType(dm_rho,dm_name);CHKERRQ(ierr);
250 // Create DM instance
251 ierr = DMMoFEMCreateMoFEM(dm_rho,&m_field,dm_name,bit_level0); CHKERRQ(ierr);
252 // Configure DM form line command options (DM itself, solvers, pre-conditioners, ... )
253 ierr = DMSetFromOptions(dm_rho); CHKERRQ(ierr);
254 // Add elements to dm (only one here)
255 ierr = DMMoFEMAddElement(dm_rho,"DENSITY_MAP_ELEMENT"); CHKERRQ(ierr);
256 // Set up problem (number dofs, partition mesh, etc.)
257 //ierr = DMSetUp(dm_rho); CHKERRQ(ierr);
258
259 {
260 ierr = DMSetUp(dm_rho); CHKERRQ(ierr);
261 // Do serial matrix
262 // const MoFEM::Problem *problem_ptr;
263 // ierr = DMMoFEMGetProblemPtr(dm_rho,&problem_ptr); CHKERRQ(ierr);
264 // ierr = m_field.build_problem(problem_ptr->getName()); CHKERRQ(ierr);
265 // // FIXME: That is not optimal partitioning
266 // ierr = m_field.partition_simple_problem(problem_ptr->getName()); CHKERRQ(ierr);
267 // ierr = m_field.partition_finite_elements(problem_ptr->getName());
268 // ierr = m_field.partition_ghost_dofs(problem_ptr->getName()); CHKERRQ(ierr);
269 }
270
271 Vec volume_vec;
272 int volume_vec_ghost[] = { 0 };
273 ierr = VecCreateGhost(
274 PETSC_COMM_WORLD,(!m_field.get_comm_rank())?1:0,1,1,volume_vec_ghost,&volume_vec
275 ); CHKERRQ(ierr);
276 ierr = VecZeroEntries(volume_vec); CHKERRQ(ierr);
277 Vec mass_vec;
278 ierr = VecDuplicate(volume_vec,&mass_vec); CHKERRQ(ierr);
279 Vec mass_vec_approx;
280 ierr = VecDuplicate(volume_vec,&mass_vec_approx); CHKERRQ(ierr);
281
282 DensityMapFe fe_volume(m_field);
283 auto &list_of_operators = fe_volume.getOpPtrVector();
284 list_of_operators.push_back(new OpVolumeCalculation("RHO",volume_vec));
285 DataFromMetaIO data_from_meta_io(meta_file);
286 ierr = data_from_meta_io.fromOptions(); CHKERRQ(ierr);
287 list_of_operators.push_back(new OpMassCalculation("RHO",mass_vec,data_from_meta_io));
288 ierr = DMoFEMLoopFiniteElements(dm_rho,"DENSITY_MAP_ELEMENT",&fe_volume); CHKERRQ(ierr);
289
290 ierr = VecAssemblyBegin(volume_vec); CHKERRQ(ierr);
291 ierr = VecAssemblyEnd(volume_vec); CHKERRQ(ierr);
292 double volume;
293 ierr = VecSum(volume_vec,&volume); CHKERRQ(ierr);
294 ierr = PetscPrintf(PETSC_COMM_WORLD,"Volume %0.32g\n",volume); CHKERRQ(ierr);
295
296 ierr = VecAssemblyBegin(mass_vec); CHKERRQ(ierr);
297 ierr = VecAssemblyEnd(mass_vec); CHKERRQ(ierr);
298 double mass;
299 ierr = VecSum(mass_vec,&mass); CHKERRQ(ierr);
300 ierr = PetscPrintf(PETSC_COMM_WORLD,"Mass %0.32g\n",mass); CHKERRQ(ierr);
301
302 ierr = VecDestroy(&volume_vec); CHKERRQ(ierr);
303 ierr = VecDestroy(&mass_vec); CHKERRQ(ierr);
304
305
306 Mat A;
307 Vec F,D;
308 ierr = DMCreateMatrix(dm_rho,&A); CHKERRQ(ierr);
309 ierr = m_field.getInterface<VecManager>()->vecCreateGhost("DM_RHO",ROW,&F); CHKERRQ(ierr);
310 ierr = m_field.getInterface<VecManager>()->vecCreateGhost("DM_RHO",COL,&D); CHKERRQ(ierr);
311
312 ierr = MatZeroEntries(A); CHKERRQ(ierr);
313 ierr = VecZeroEntries(F); CHKERRQ(ierr);
314
315 CommonData common_data;
316 common_data.globA = A;
317 common_data.globF = F;
318
319 common_data.distanceValue = boost::shared_ptr<VectorDouble>(new VectorDouble());
320 common_data.distanceGradient = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
321 boost::shared_ptr<VectorDouble> dist_value = common_data.distanceValue;
322 boost::shared_ptr<MatrixDouble> grad_value = common_data.distanceGradient;
323
324 DensityMapFe fe_density_approximation(m_field);
325
326 fe_density_approximation.getOpPtrVector().push_back
327 (new OpCalculateScalarFieldValues("DISTANCE_FROM_SURFACE",dist_value));
328 fe_density_approximation.getOpPtrVector().push_back
329 (new OpCalculateScalarFieldGradient<3>("DISTANCE_FROM_SURFACE",grad_value));
330 fe_density_approximation.getOpPtrVector().push_back
331 (new OpCalculateLhs("RHO","RHO",common_data,data_from_meta_io));
332 fe_density_approximation.getOpPtrVector().push_back
333 (new OpCalulatefRhoAtGaussPts("RHO",common_data,data_from_meta_io,surface_distance));
334 fe_density_approximation.getOpPtrVector().push_back
335 (new OpAssmbleRhs("RHO",common_data));
336 ierr = DMoFEMLoopFiniteElements(dm_rho,"DENSITY_MAP_ELEMENT",&fe_density_approximation); CHKERRQ(ierr);
337
338 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
339 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
340 ierr = VecGhostUpdateBegin(F,ADD_VALUES,SCATTER_REVERSE);
341 ierr = VecGhostUpdateEnd(F,ADD_VALUES,SCATTER_REVERSE);
342 ierr = VecAssemblyBegin(F); CHKERRQ(ierr);
343 ierr = VecAssemblyEnd(F); CHKERRQ(ierr);
344
345 // MatView(A,PETSC_VIEWER_DRAW_WORLD);
346 // std::string wait;
347 // std::cin >> wait;
348
349 {
350 KSP ksp;
351 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp); CHKERRQ(ierr);
352 ierr = KSPSetOperators(ksp,A,A); CHKERRQ(ierr);
353 ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);
354 ierr = KSPSolve(ksp,F,D); CHKERRQ(ierr); CHKERRQ(ierr);
355 }
356
357 ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
358 ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
359 ierr = m_field.getInterface<VecManager>()->setGlobalGhostVector("DM_RHO",COL,D,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
360
361 DensityMapFe fe_post_process_calulations(m_field);
362
363 fe_post_process_calulations.getOpPtrVector().push_back(
364 new OpMassCalculationFromApprox("RHO",mass_vec_approx)
365 );
366 ierr = DMoFEMLoopFiniteElements(dm_rho,"DENSITY_MAP_ELEMENT",&fe_post_process_calulations); CHKERRQ(ierr);
367
368 ierr = VecAssemblyBegin(mass_vec_approx); CHKERRQ(ierr);
369 ierr = VecAssemblyEnd(mass_vec_approx); CHKERRQ(ierr);
370 double mass_approx;
371 ierr = VecSum(mass_vec_approx,&mass_approx); CHKERRQ(ierr);
372 ierr = PetscPrintf(PETSC_COMM_WORLD,"Mass after approximation %0.32g\n",mass_approx); CHKERRQ(ierr);
373 if(is_atom_test) {
374 if( abs(mass_approx - volume) > 1.0e-2) {
375 SETERRQ(PETSC_COMM_SELF,MOFEM_ATOM_TEST_INVALID,"atom test diverged");
376 }
377 }
378 ierr = VecDestroy(&mass_vec_approx); CHKERRQ(ierr);
379
380 PostProcVolumeOnRefinedMesh post_proc(m_field);
381 ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
382 ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
383 ierr = post_proc.addFieldValuesPostProc("RHO"); CHKERRQ(ierr);
384 ierr = post_proc.addFieldValuesPostProc("DISTANCE_FROM_SURFACE"); CHKERRQ(ierr);
385 ierr = post_proc.addFieldValuesGradientPostProc("RHO"); CHKERRQ(ierr);
386 ierr = post_proc.addFieldValuesGradientPostProc("DISTANCE_FROM_SURFACE"); CHKERRQ(ierr);
387 ierr = DMoFEMLoopFiniteElements(dm_rho,"DENSITY_MAP_ELEMENT",&post_proc); CHKERRQ(ierr);
388 ierr = post_proc.writeFile("out.h5m"); CHKERRQ(ierr);
389
390 ierr = VecDestroy(&F); CHKERRQ(ierr);
391 ierr = VecDestroy(&D); CHKERRQ(ierr);
392 ierr = MatDestroy(&A); CHKERRQ(ierr);
393 ierr = DMDestroy(&dm_rho);
394
395 // Delete element which will be no longer used
396 ierr = m_field.delete_finite_element("DENSITY_MAP_ELEMENT",2); CHKERRQ(ierr);
397 if (m_field.get_comm_rank() == 0)
398 CHKERR moab.write_file("analysis_mesh.h5m");
399
400
401 } CATCH_ERRORS;
402
404
405 return (0);
406}
Operator can be used with any volume element to calculate sum of volumes of all volumes in the set.
Post-process fields on refined mesh.
int main()
Definition: adol-c_atom.cpp:46
cholesky decomposition
@ COL
Definition: definitions.h:123
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ H1
continuous field
Definition: definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
@ BLOCKSET
Definition: definitions.h:148
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:454
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
#define CHKERR
Inline error check.
Definition: definitions.h:535
@ 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 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
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.
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
MoFEMErrorCode writeFile(const std::string file_name, const char *file_type="MOAB", const char *file_options="PARALLEL=WRITE_PART")
wrote results in (MOAB) format, use "file_name.h5m"
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.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
double D
static char help[]
char mesh_file_name[255]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
constexpr AssemblyType A
Load data from MetaImage file, translate grayscale values to densities.
Assemble RHS vector f.
Assemble LHS matrix K.
Assemble local vector containing density data.
Calculate mass after approximation.
Calculate mass before approximation.
Create KDTree on surface of the mesh and calculate distance.
Definition: DensityMaps.hpp:32
PetscErrorCode setDistanceFromSurface(const std::string field_name)
Set distance to vertices on mesh.
PetscErrorCode takeASkin(Range &volume)
Take a skin from a mesh.
Definition: DensityMaps.hpp:55
PetscErrorCode buildTree()
Build tree.
Definition: DensityMaps.hpp:72
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode delete_finite_element(const std::string name, int verb=DEFAULT_VERBOSITY)=0
delete finite element from mofem database
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.
virtual int get_comm_rank() const =0
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.
DEPRECATED MoFEMErrorCode seed_ref_level_3D(const EntityHandle meshset, const BitRefLevel &bit, int verb=-1)
seed 2D entities in the meshset and their adjacencies (only TETs adjacencies) in a particular BitRefL...
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Interface for managing meshsets containing materials and boundary conditions.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Projection of edge entities with one mid-node on hierarchical basis.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Post processing.