v0.10.0
Functions | Variables
map_disp_prism.cpp File Reference
#include <MoFEM.hpp>
#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#include <algorithm>
#include <iterator>
#include <boost/tokenizer.hpp>
#include <PostProcOnRefMesh.hpp>
#include <DispMap.hpp>

Go to the source code of this file.

Functions

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

Variables

static char help [] = "\n"
 

Function Documentation

◆ main()

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

Definition at line 35 of file map_disp_prism.cpp.

35  {
36  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
37 
38  try {
39 
40  moab::Core mb_instance;
41  moab::Interface &moab = mb_instance;
42 
43  // global variables
44  char mesh_file_name[255];
45  char data_file_name1[255];
46  char data_file_name2[255];
47  PetscBool flg_file = PETSC_FALSE;
48  PetscInt order = 1;
49  PetscBool flg_n_part = PETSC_FALSE;
50  PetscInt n_partas = 1;
51  PetscBool flg_thickness = PETSC_FALSE;
52  double thickness = 0;
53 
54  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "none", "none");
55  CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
56  mesh_file_name, 255, &flg_file);
57 
58  CHKERR PetscOptionsString("-my_data_x", "data file name", "", "data_x.h5m",
59  data_file_name1, 255, PETSC_NULL);
60 
61  CHKERR PetscOptionsString("-my_data_y", "data file name", "", "data_y.h5m",
62  data_file_name2, 255, PETSC_NULL);
63 
64  CHKERR PetscOptionsInt("-my_order", "default approximation order", "", 1,
65  &order, PETSC_NULL);
66 
67  CHKERR PetscOptionsInt("-my_nparts", "number of parts", "", 1, &n_partas,
68  &flg_n_part);
69 
70  CHKERR PetscOptionsScalar("-my_thinckness", "top layer thickness", "",
71  thickness, &thickness, &flg_thickness);
72 
73  ierr = PetscOptionsEnd();
74  CHKERRQ(ierr);
75 
76  if (flg_file != PETSC_TRUE) {
77  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
78  }
79 
80  if (flg_n_part != PETSC_TRUE) {
81  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR partitioning number not given");
82  }
83 
84  if (flg_thickness != PETSC_TRUE) {
85  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR thickness of top layer not given");
86  }
87 
88  const char *option;
89  option = "";
90  CHKERR moab.load_file(mesh_file_name, 0, option);
91  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
92  if (pcomm == NULL)
93  pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
94 
95  MoFEM::Core core(moab);
96  MoFEM::Interface &m_field = core;
97  // meshset consisting all entities in mesh
98  EntityHandle root_set = moab.get_root_set();
99 
100  // Seed all mesh entities to MoFEM database, those entities can be
101  // potentially used as finite elements or as entities which carry some
102  // approximation field.
103  BitRefLevel bit_level0;
104  bit_level0.set(0);
105  CHKERR m_field.seed_ref_level_3D(root_set, bit_level0);
106  // Add field
107  CHKERR m_field.add_field("DISP_X", H1, AINSWORTH_LEGENDRE_BASE, 1);
108  CHKERR m_field.add_field("DISP_Y", H1, AINSWORTH_LEGENDRE_BASE, 1);
109 
110  MeshsetsManager *meshsets_manager_ptr;
111  CHKERR m_field.query_interface(meshsets_manager_ptr);
112 
113  Range surface;
114  CHKERR meshsets_manager_ptr->getEntitiesByDimension(101, SIDESET, 2,
115  surface, true);
116  CHKERR m_field.add_ents_to_field_by_type(surface, MBTRI, "DISP_X");
117  CHKERR m_field.add_ents_to_field_by_type(surface, MBTRI, "DISP_Y");
118 
119  CHKERR m_field.set_field_order(0, MBVERTEX, "DISP_X", 1);
120  CHKERR m_field.set_field_order(0, MBEDGE, "DISP_X", order);
121  CHKERR m_field.set_field_order(0, MBTRI, "DISP_X", order);
122 
123  CHKERR m_field.set_field_order(0, MBVERTEX, "DISP_Y", 1);
124  CHKERR m_field.set_field_order(0, MBEDGE, "DISP_Y", order);
125  CHKERR m_field.set_field_order(0, MBTRI, "DISP_Y", order);
126 
127  CHKERR m_field.build_fields();
128 
129  CHKERR m_field.add_finite_element("DISP_MAP_ELEMENT");
130  CHKERR m_field.modify_finite_element_add_field_row("DISP_MAP_ELEMENT",
131  "DISP_X");
132  CHKERR m_field.modify_finite_element_add_field_col("DISP_MAP_ELEMENT",
133  "DISP_X");
134  CHKERR m_field.modify_finite_element_add_field_data("DISP_MAP_ELEMENT",
135  "DISP_X");
136  CHKERR m_field.modify_finite_element_add_field_data("DISP_MAP_ELEMENT",
137  "DISP_Y");
138 
139  // Add entities to that element
140  CHKERR m_field.add_ents_to_finite_element_by_type(surface, MBTRI,
141  "DISP_MAP_ELEMENT");
142 
143  // build finite elements
144  CHKERR m_field.build_finite_elements();
145  // build adjacencies between elements and degrees of freedom
146  CHKERR m_field.build_adjacencies(bit_level0);
147 
148  // set discrete manager
149  DM dm_disp;
150  DMType dm_name = "DM_DISP";
151  // Register DM problem
152  CHKERR DMRegister_MoFEM(dm_name);
153  CHKERR DMCreate(PETSC_COMM_WORLD, &dm_disp);
154  CHKERR DMSetType(dm_disp, dm_name);
155  // Create DM instance
156  CHKERR DMMoFEMCreateMoFEM(dm_disp, &m_field, dm_name, bit_level0);
157  // Configure DM form line command options (DM itself, solvers,
158  // pre-conditioners, ... )
159  CHKERR DMSetFromOptions(dm_disp);
160  // Add elements to dm (only one here)
161  CHKERR DMMoFEMAddElement(dm_disp, "DISP_MAP_ELEMENT");
162  // Set up problem (number dofs, partition mesh, etc.)
163  CHKERR DMSetUp(dm_disp);
164 
165  DispMapFe fe_face(m_field);
166 
167  // string file(argv[1]); cout << "argument 1:" << argv[1] <<endl;
168  // string file2(argv[2]); cout << "argument 2:" << argv[2] <<endl;
169 
170  string file1(data_file_name1);
171  DataFromFiles data_from_files(file1);
172  CHKERR data_from_files.fromOptions();
173  CHKERR data_from_files.loadFileData();
174 
175  string file2(data_file_name2);
176  DataFromFiles data_from_files2(file2);
177  CHKERR data_from_files2.fromOptions();
178  CHKERR data_from_files2.loadFileData();
179 
180  boost::ptr_vector<MoFEM::ForcesAndSourcesCore::UserDataOperator>
181  &list_of_operators = fe_face.getOpPtrVector();
182 
183  Mat A;
184  Vec Fx, Dx;
185  CHKERR DMCreateMatrix(dm_disp, &A);
186  CHKERR m_field.query_interface<VecManager>()->vecCreateGhost("DM_DISP", ROW,
187  &Fx);
188  CHKERR m_field.query_interface<VecManager>()->vecCreateGhost("DM_DISP", COL,
189  &Dx);
190  Vec Fy, Dy;
191  CHKERR VecDuplicate(Fx, &Fy);
192  CHKERR VecDuplicate(Dx, &Dy);
193 
194  CHKERR MatZeroEntries(A);
195  CHKERR VecZeroEntries(Fx);
196  CHKERR VecZeroEntries(Fy);
197 
198  CommonData common_data;
199  common_data.globA = A;
200 
201  common_data.globF = Fx;
202  list_of_operators.push_back(
203  new OpCalculateInvJacForFace(common_data.invJac));
204  list_of_operators.push_back(new OpSetInvJacH1ForFace(common_data.invJac));
205  list_of_operators.push_back(
206  new OpCalculateLhs("DISP_X", "DISP_X", common_data, data_from_files));
207  list_of_operators.push_back(
208  new OpCalulatefRhoAtGaussPts("DISP_X", common_data, data_from_files));
209  list_of_operators.push_back(new OpAssmbleRhs("DISP_X", common_data));
210 
211  CHKERR DMoFEMLoopFiniteElements(dm_disp, "DISP_MAP_ELEMENT", &fe_face);
212 
213  // cout << "error.no !" <<error++ << endl;
214  CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
215  CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
216  CHKERR VecGhostUpdateBegin(Fx, ADD_VALUES, SCATTER_REVERSE);
217  CHKERR VecGhostUpdateEnd(Fx, ADD_VALUES, SCATTER_REVERSE);
218  CHKERR VecAssemblyBegin(Fx);
219  CHKERR VecAssemblyEnd(Fx);
220 
221  // int error=1;
222  common_data.globF = Fy;
223  list_of_operators.clear();
224  list_of_operators.push_back(
225  new OpCalculateInvJacForFace(common_data.invJac));
226  list_of_operators.push_back(new OpSetInvJacH1ForFace(common_data.invJac));
227  list_of_operators.push_back(
228  new OpCalulatefRhoAtGaussPts("DISP_X", common_data, data_from_files2));
229  list_of_operators.push_back(new OpAssmbleRhs("DISP_X", common_data));
230  // cout << "error.no !"<< error++<< endl;
231 
232  CHKERR DMoFEMLoopFiniteElements(dm_disp, "DISP_MAP_ELEMENT", &fe_face);
233  // cout << "error.no !" <<error++ << endl;
234  CHKERR VecGhostUpdateBegin(Fy, ADD_VALUES, SCATTER_REVERSE);
235  CHKERR VecGhostUpdateEnd(Fy, ADD_VALUES, SCATTER_REVERSE);
236  CHKERR VecAssemblyBegin(Fy);
237  CHKERR VecAssemblyEnd(Fy);
238 
239  // VecView(F,PETSC_VIEWER_DRAW_WORLD);
240  // MatView(A,PETSC_VIEWER_DRAW_WORLD);
241  // std::string wait;
242  // std::cin >> wait;
243 
244  // solve
245  {
246  KSP ksp;
247  CHKERR KSPCreate(PETSC_COMM_WORLD, &ksp);
248  CHKERR KSPSetOperators(ksp, A, A);
249  CHKERR KSPSetFromOptions(ksp);
250  CHKERR KSPSolve(ksp, Fx, Dx);
251  CHKERR KSPSolve(ksp, Fy, Dy);
252  }
253 
254  CHKERR VecGhostUpdateBegin(Dx, INSERT_VALUES, SCATTER_FORWARD);
255  CHKERR VecGhostUpdateEnd(Dx, INSERT_VALUES, SCATTER_FORWARD);
256  CHKERR m_field.query_interface<VecManager>()->setGlobalGhostVector(
257  "DM_DISP", COL, Dx, INSERT_VALUES, SCATTER_REVERSE);
258  CHKERR VecGhostUpdateBegin(Dy, INSERT_VALUES, SCATTER_FORWARD);
259  CHKERR VecGhostUpdateEnd(Dy, INSERT_VALUES, SCATTER_FORWARD);
260  CHKERR m_field.query_interface<VecManager>()->setOtherGlobalGhostVector(
261  "DM_DISP", "DISP_X", "DISP_Y", COL, Dy, INSERT_VALUES, SCATTER_REVERSE);
262 
263  // Remove rigid translation
264  double sum_dx = 0;
265  int nb_dx = 0;
266  for (_IT_GET_DOFS_FIELD_BY_NAME_AND_TYPE_FOR_LOOP_(m_field, "DISP_X",
267  MBVERTEX, dof)) {
268  sum_dx += (*dof)->getFieldData();
269  nb_dx++;
270  }
271  sum_dx /= nb_dx;
272  for (_IT_GET_DOFS_FIELD_BY_NAME_AND_TYPE_FOR_LOOP_(m_field, "DISP_X",
273  MBVERTEX, dof)) {
274  (*dof)->getFieldData() -= sum_dx;
275  }
276  double sum_dy = 0;
277  int nb_dy = 0;
278  for (_IT_GET_DOFS_FIELD_BY_NAME_AND_TYPE_FOR_LOOP_(m_field, "DISP_Y",
279  MBVERTEX, dof)) {
280  sum_dy += (*dof)->getFieldData();
281  nb_dy++;
282  }
283  sum_dy /= nb_dy;
284  for (_IT_GET_DOFS_FIELD_BY_NAME_AND_TYPE_FOR_LOOP_(m_field, "DISP_Y",
285  MBVERTEX, dof)) {
286  (*dof)->getFieldData() -= sum_dy;
287  }
288 
290  CHKERR post_proc.generateReferenceElementMesh();
291  CHKERR post_proc.addFieldValuesPostProc("DISP_X");
292  CHKERR post_proc.addFieldValuesPostProc("DISP_Y");
293  CHKERR DMoFEMLoopFiniteElements(dm_disp, "DISP_MAP_ELEMENT", &post_proc);
294 
295  CHKERR post_proc.postProcMesh.write_file("out_disp.h5m", "MOAB",
296  "PARALLEL=WRITE_PART");
297 
298  CHKERR VecDestroy(&Fx);
299  CHKERR VecDestroy(&Fy);
300  CHKERR VecDestroy(&Dx);
301  CHKERR VecDestroy(&Dy);
302  CHKERR MatDestroy(&A);
303  CHKERR DMDestroy(&dm_disp);
304 
305  // Delete element which will be no longer used
306  CHKERR m_field.delete_finite_element("DISP_MAP_ELEMENT");
307 
308  // Add top layer
309  {
310  MeshsetsManager *mmanager_ptr;
311  CHKERR m_field.query_interface(mmanager_ptr);
312  Range ents_2nd_layer;
313  CHKERR mmanager_ptr->getEntitiesByDimension(101, SIDESET, 2,
314  ents_2nd_layer, true);
315  Range tris_nodes;
316  CHKERR moab.get_connectivity(ents_2nd_layer, tris_nodes, true);
317  std::map<EntityHandle, EntityHandle> map_nodes;
318  for (Range::iterator nit = tris_nodes.begin(); nit != tris_nodes.end();
319  nit++) {
320  double coords[3];
321  CHKERR moab.get_coords(&*nit, 1, coords);
322  coords[2] += thickness;
323  CHKERR moab.create_vertex(coords, map_nodes[*nit]);
324  }
325 
326  Range prisms;
327  Range top_tris;
328  for (Range::iterator tit = ents_2nd_layer.begin();
329  tit != ents_2nd_layer.end(); tit++) {
330  int num_nodes;
331  const EntityHandle *conn;
332  CHKERR moab.get_connectivity(*tit, conn, num_nodes);
333  EntityHandle conn_prism[6];
334  for (int n = 0; n != 3; n++) {
335  conn_prism[n] = conn[n];
336  conn_prism[n + 3] = map_nodes[conn[n]];
337  }
338  EntityHandle prism;
339  CHKERR moab.create_element(MBPRISM, conn_prism, 6, prism);
340  Range adj;
341  CHKERR moab.get_adjacencies(&prism, 1, 2, true, adj);
342  CHKERR moab.get_adjacencies(&prism, 1, 1, true, adj);
343  prisms.insert(prism);
344  EntityHandle tri;
345  CHKERR moab.side_element(prism, 2, 4, tri);
346  top_tris.insert(tri);
347  }
348  CHKERR mmanager_ptr->addMeshset(SIDESET, 202);
349  CHKERR mmanager_ptr->addEntitiesToMeshset(SIDESET, 202, top_tris);
350  CHKERR mmanager_ptr->addMeshset(BLOCKSET, 2, "MAT_ELASTIC_2");
351  CHKERR mmanager_ptr->addEntitiesToMeshset(BLOCKSET, 2, prisms);
352  Mat_Elastic mat_elastic;
353  mat_elastic.data.Young = 1;
354  mat_elastic.data.Poisson = 0;
356  mat_elastic);
357  }
358 
359  {
360  // Add one node to fix
361  MeshsetsManager *mmanager_ptr;
362  CHKERR m_field.query_interface(mmanager_ptr);
363  Range ents_1st_layer;
364  CHKERR mmanager_ptr->getEntitiesByDimension(202, SIDESET, 2,
365  ents_1st_layer, true);
366  int num_nodes;
367  const EntityHandle *conn;
368  CHKERR moab.get_connectivity(ents_1st_layer[0], conn, num_nodes);
369  Skinner skin(&moab);
370  Range skin_edges;
371  CHKERR skin.find_skin(0, ents_1st_layer, false, skin_edges);
372  Range edges;
373  CHKERR moab.get_adjacencies(&*ents_1st_layer.begin(), 1, 1, false, edges);
374  EntityHandle meshset;
375  CHKERR mmanager_ptr->getMeshset(202, SIDESET, meshset);
376  CHKERR moab.add_entities(meshset, conn, 1);
377  CHKERR moab.add_entities(meshset, skin_edges);
378  }
379 
380  {
381  Range ents3d;
382  CHKERR moab.get_entities_by_dimension(0, 3, ents3d, false);
383  CHKERR m_field.partition_mesh(ents3d, 3, 2, n_partas);
384  }
385 
386  CHKERR moab.write_file("analysis_mesh.h5m");
387  }
388  CATCH_ERRORS;
389 
391 
392  PetscFunctionReturn(0);
393 }
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
MoFEMErrorCode getMeshset(const int ms_id, const unsigned int cubit_bc_type, EntityHandle &meshset) const
get meshset from CUBIT Id and CUBIT type
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
Deprecated interface functions.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
DEPRECATED MoFEMErrorCode partition_mesh(const Range &ents, const int dim, const int adj_dim, const int n_parts, int verb=-1)
Set partition tag to each finite element in the problem.
CoreTmp< 0 > Core
Definition: Core.hpp:1129
#define _IT_GET_DOFS_FIELD_BY_NAME_AND_TYPE_FOR_LOOP_(MFIELD, NAME, TYPE, IT)
loop over all dofs from a moFEM field and particular field
Definition: Interface.hpp:1893
auto post_proc
MoFEMErrorCode getEntitiesByDimension(const int ms_id, const unsigned int cubit_bc_type, const int dimension, Range &entities, const bool recursive=true) const
get entities from CUBIT/meshset of a particular entity dimensionNodeset can contain nodes,...
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:425
DEPRECATED MoFEMErrorCode query_interface(IFace *&ptr) const
Elastic material data structure.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:60
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.
Interface for managing meshsets containing materials and boundary conditions.
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.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
char mesh_file_name[255]
static Index< 'n', 3 > n
static char help[]
constexpr int order
Calculate inverse of jacobian for face element.
MoFEMErrorCode addMeshset(const CubitBCType cubit_bc_type, const int ms_id, const std::string name="")
add cubit meshset
CHKERRQ(ierr)
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.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
MoFEMErrorCode addEntitiesToMeshset(const CubitBCType cubit_bc_type, const int ms_id, const Range &ents)
add entities to cubit meshset
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:36
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: DMMMoFEM.cpp:105
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
#define CHKERR
Inline error check.
Definition: definitions.h:604
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
Transform local reference derivatives of shape functions to global derivatives.
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:292
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
Postprocess on face.
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
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
Core (interface) class.
Definition: Core.hpp:77
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:441
continuous field
Definition: definitions.h:177
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...
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
MoFEMErrorCode setAtributesByDataStructure(const CubitBCType cubit_bc_type, const int ms_id, const GenericAttributeData &data, const std::string name="")
set (material) data structure to cubit meshset
Mat globA
Global matrix.
Definition: DispMap.hpp:179
virtual MoFEMErrorCode delete_finite_element(const std::string name, int verb=DEFAULT_VERBOSITY)=0
delete finite element from mofem database

Variable Documentation

◆ help

char help[] = "\n"
static

Definition at line 33 of file map_disp_prism.cpp.