v0.10.0
Functions | Variables
map_disp.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.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 
52  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "none", "none");
53  CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
54  mesh_file_name, 255, &flg_file);
55 
56  CHKERR PetscOptionsString("-my_data_x", "data file name", "", "data_x.h5m",
57  data_file_name1, 255, PETSC_NULL);
58 
59  CHKERR PetscOptionsString("-my_data_y", "data file name", "", "data_y.h5m",
60  data_file_name2, 255, PETSC_NULL);
61 
62  CHKERR PetscOptionsInt("-my_order", "default approximation order", "", 1,
63  &order, PETSC_NULL);
64 
65  CHKERR PetscOptionsInt("-my_nparts", "number of parts", "", 1, &n_partas,
66  &flg_n_part);
67 
68  ierr = PetscOptionsEnd();
69  CHKERRQ(ierr);
70 
71  if (flg_file != PETSC_TRUE) {
72  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
73  }
74 
75  if (flg_n_part != PETSC_TRUE) {
76  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR partitioning number not given");
77  }
78 
79  const char *option;
80  option = "";
81  CHKERR moab.load_file(mesh_file_name, 0, option);
82  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
83  if (pcomm == NULL)
84  pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
85 
86  MoFEM::Core core(moab);
87  MoFEM::Interface &m_field = core;
88  // meshset consisting all entities in mesh
89  EntityHandle root_set = moab.get_root_set();
90 
91  // Seed all mesh entities to MoFEM database, those entities can be
92  // potentially used as finite elements or as entities which carry some
93  // approximation field.
94  BitRefLevel bit_level0;
95  bit_level0.set(0);
96  CHKERR m_field.seed_ref_level_3D(root_set, bit_level0);
97  // Add field
98  CHKERR m_field.add_field("DISP_X", H1, AINSWORTH_LEGENDRE_BASE, 1);
99  CHKERR m_field.add_field("DISP_Y", H1, AINSWORTH_LEGENDRE_BASE, 1);
100 
101  MeshsetsManager *meshsets_manager_ptr;
102  CHKERR m_field.query_interface(meshsets_manager_ptr);
103 
104  Range surface;
105  CHKERR meshsets_manager_ptr->getEntitiesByDimension(101, SIDESET, 2,
106  surface, true);
107  CHKERR m_field.add_ents_to_field_by_type(surface, MBTRI, "DISP_X");
108  CHKERR m_field.add_ents_to_field_by_type(surface, MBTRI, "DISP_Y");
109 
110  CHKERR m_field.set_field_order(0, MBVERTEX, "DISP_X", 1);
111  CHKERR m_field.set_field_order(0, MBEDGE, "DISP_X", order);
112  CHKERR m_field.set_field_order(0, MBTRI, "DISP_X", order);
113 
114  CHKERR m_field.set_field_order(0, MBVERTEX, "DISP_Y", 1);
115  CHKERR m_field.set_field_order(0, MBEDGE, "DISP_Y", order);
116  CHKERR m_field.set_field_order(0, MBTRI, "DISP_Y", order);
117 
118  CHKERR m_field.build_fields();
119 
120  CHKERR m_field.add_finite_element("DISP_MAP_ELEMENT");
121  CHKERR m_field.modify_finite_element_add_field_row("DISP_MAP_ELEMENT",
122  "DISP_X");
123  CHKERR m_field.modify_finite_element_add_field_col("DISP_MAP_ELEMENT",
124  "DISP_X");
125  CHKERR m_field.modify_finite_element_add_field_data("DISP_MAP_ELEMENT",
126  "DISP_X");
127  CHKERR m_field.modify_finite_element_add_field_data("DISP_MAP_ELEMENT",
128  "DISP_Y");
129 
130  // Add entities to that element
131  CHKERR m_field.add_ents_to_finite_element_by_type(surface, MBTRI,
132  "DISP_MAP_ELEMENT");
133 
134  // build finite elements
135  CHKERR m_field.build_finite_elements();
136  // build adjacencies between elements and degrees of freedom
137  CHKERR m_field.build_adjacencies(bit_level0);
138 
139  // set discrete manager
140  DM dm_disp;
141  DMType dm_name = "DM_DISP";
142  // Register DM problem
143  CHKERR DMRegister_MoFEM(dm_name);
144  CHKERR DMCreate(PETSC_COMM_WORLD, &dm_disp);
145  CHKERR DMSetType(dm_disp, dm_name);
146  // Create DM instance
147  CHKERR DMMoFEMCreateMoFEM(dm_disp, &m_field, dm_name, bit_level0);
148  // Configure DM form line command options (DM itself, solvers,
149  // pre-conditioners, ... )
150  CHKERR DMSetFromOptions(dm_disp);
151  // Add elements to dm (only one here)
152  CHKERR DMMoFEMAddElement(dm_disp, "DISP_MAP_ELEMENT");
153  // Set up problem (number dofs, partition mesh, etc.)
154  CHKERR DMSetUp(dm_disp);
155 
156  DispMapFe fe_face(m_field);
157 
158  // string file(argv[1]); cout << "argument 1:" << argv[1] <<endl;
159  // string file2(argv[2]); cout << "argument 2:" << argv[2] <<endl;
160 
161  string file1(data_file_name1);
162  DataFromFiles data_from_files(file1);
163  CHKERR data_from_files.fromOptions();
164  CHKERR data_from_files.loadFileData();
165 
166  string file2(data_file_name2);
167  DataFromFiles data_from_files2(file2);
168  CHKERR data_from_files2.fromOptions();
169  CHKERR data_from_files2.loadFileData();
170 
171  boost::ptr_vector<MoFEM::ForcesAndSourcesCore::UserDataOperator>
172  &list_of_operators = fe_face.getOpPtrVector();
173 
174  Mat A;
175  Vec Fx, Dx;
176  CHKERR DMCreateMatrix(dm_disp, &A);
177  CHKERR m_field.query_interface<VecManager>()->vecCreateGhost("DM_DISP", ROW,
178  &Fx);
179  CHKERR m_field.query_interface<VecManager>()->vecCreateGhost("DM_DISP", COL,
180  &Dx);
181  Vec Fy, Dy;
182  CHKERR VecDuplicate(Fx, &Fy);
183  CHKERR VecDuplicate(Dx, &Dy);
184 
185  CHKERR MatZeroEntries(A);
186  CHKERR VecZeroEntries(Fx);
187  CHKERR VecZeroEntries(Fy);
188 
189  CommonData common_data;
190  common_data.globA = A;
191 
192  common_data.globF = Fx;
193  list_of_operators.push_back(
194  new OpCalculateInvJacForFace(common_data.invJac));
195  list_of_operators.push_back(new OpSetInvJacH1ForFace(common_data.invJac));
196  list_of_operators.push_back(
197  new OpCalculateLhs("DISP_X", "DISP_X", common_data, data_from_files));
198  list_of_operators.push_back(
199  new OpCalulatefRhoAtGaussPts("DISP_X", common_data, data_from_files));
200  list_of_operators.push_back(new OpAssmbleRhs("DISP_X", common_data));
201 
202  CHKERR DMoFEMLoopFiniteElements(dm_disp, "DISP_MAP_ELEMENT", &fe_face);
203 
204  // cout << "error.no !" <<error++ << endl;
205  CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
206  CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
207  CHKERR VecGhostUpdateBegin(Fx, ADD_VALUES, SCATTER_REVERSE);
208  CHKERR VecGhostUpdateEnd(Fx, ADD_VALUES, SCATTER_REVERSE);
209  CHKERR VecAssemblyBegin(Fx);
210  CHKERR VecAssemblyEnd(Fx);
211 
212  // int error=1;
213  common_data.globF = Fy;
214  list_of_operators.clear();
215  list_of_operators.push_back(
216  new OpCalculateInvJacForFace(common_data.invJac));
217  list_of_operators.push_back(new OpSetInvJacH1ForFace(common_data.invJac));
218  list_of_operators.push_back(
219  new OpCalulatefRhoAtGaussPts("DISP_X", common_data, data_from_files2));
220  list_of_operators.push_back(new OpAssmbleRhs("DISP_X", common_data));
221  // cout << "error.no !"<< error++<< endl;
222 
223  CHKERR DMoFEMLoopFiniteElements(dm_disp, "DISP_MAP_ELEMENT", &fe_face);
224  // cout << "error.no !" <<error++ << endl;
225  CHKERR VecGhostUpdateBegin(Fy, ADD_VALUES, SCATTER_REVERSE);
226  CHKERR VecGhostUpdateEnd(Fy, ADD_VALUES, SCATTER_REVERSE);
227  CHKERR VecAssemblyBegin(Fy);
228  CHKERR VecAssemblyEnd(Fy);
229 
230  // VecView(F,PETSC_VIEWER_DRAW_WORLD);
231  // MatView(A,PETSC_VIEWER_DRAW_WORLD);
232  // std::string wait;
233  // std::cin >> wait;
234 
235  // solve
236  {
237  KSP ksp;
238  CHKERR KSPCreate(PETSC_COMM_WORLD, &ksp);
239  CHKERR KSPSetOperators(ksp, A, A);
240  CHKERR KSPSetFromOptions(ksp);
241  CHKERR KSPSolve(ksp, Fx, Dx);
242  CHKERR KSPSolve(ksp, Fy, Dy);
243  }
244 
245  CHKERR VecGhostUpdateBegin(Dx, INSERT_VALUES, SCATTER_FORWARD);
246  CHKERR VecGhostUpdateEnd(Dx, INSERT_VALUES, SCATTER_FORWARD);
247  CHKERR m_field.query_interface<VecManager>()->setGlobalGhostVector(
248  "DM_DISP", COL, Dx, INSERT_VALUES, SCATTER_REVERSE);
249  CHKERR VecGhostUpdateBegin(Dy, INSERT_VALUES, SCATTER_FORWARD);
250  CHKERR VecGhostUpdateEnd(Dy, INSERT_VALUES, SCATTER_FORWARD);
251  CHKERR m_field.query_interface<VecManager>()->setOtherGlobalGhostVector(
252  "DM_DISP", "DISP_X", "DISP_Y", COL, Dy, INSERT_VALUES, SCATTER_REVERSE);
253 
254  // Remove rigid translation
255  double sum_dx = 0;
256  int nb_dx = 0;
257  for (_IT_GET_DOFS_FIELD_BY_NAME_AND_TYPE_FOR_LOOP_(m_field, "DISP_X",
258  MBVERTEX, dof)) {
259  sum_dx += (*dof)->getFieldData();
260  nb_dx++;
261  }
262  sum_dx /= nb_dx;
263  for (_IT_GET_DOFS_FIELD_BY_NAME_AND_TYPE_FOR_LOOP_(m_field, "DISP_X",
264  MBVERTEX, dof)) {
265  (*dof)->getFieldData() -= sum_dx;
266  }
267  double sum_dy = 0;
268  int nb_dy = 0;
269  for (_IT_GET_DOFS_FIELD_BY_NAME_AND_TYPE_FOR_LOOP_(m_field, "DISP_Y",
270  MBVERTEX, dof)) {
271  sum_dy += (*dof)->getFieldData();
272  nb_dy++;
273  }
274  sum_dy /= nb_dy;
275  for (_IT_GET_DOFS_FIELD_BY_NAME_AND_TYPE_FOR_LOOP_(m_field, "DISP_Y",
276  MBVERTEX, dof)) {
277  (*dof)->getFieldData() -= sum_dy;
278  }
279 
281  CHKERR post_proc.generateReferenceElementMesh();
282  CHKERR post_proc.addFieldValuesPostProc("DISP_X");
283  CHKERR post_proc.addFieldValuesPostProc("DISP_Y");
284  CHKERR DMoFEMLoopFiniteElements(dm_disp, "DISP_MAP_ELEMENT", &post_proc);
285 
286  CHKERR post_proc.postProcMesh.write_file("out_disp.h5m", "MOAB",
287  "PARALLEL=WRITE_PART");
288 
289  CHKERR VecDestroy(&Fx);
290  CHKERR VecDestroy(&Fy);
291  CHKERR VecDestroy(&Dx);
292  CHKERR VecDestroy(&Dy);
293  CHKERR MatDestroy(&A);
294  CHKERR DMDestroy(&dm_disp);
295 
296  // Delete element which will be no longer used
297  CHKERR m_field.delete_finite_element("DISP_MAP_ELEMENT");
298 
299  {
300  // Add one node to fix
301  MeshsetsManager *mmanager_ptr;
302  CHKERR m_field.query_interface(mmanager_ptr);
303  Range ents_1st_layer;
304  CHKERR mmanager_ptr->getEntitiesByDimension(202, SIDESET, 2,
305  ents_1st_layer, true);
306  int num_nodes;
307  const EntityHandle *conn;
308  CHKERR moab.get_connectivity(ents_1st_layer[0], conn, num_nodes);
309  Skinner skin(&moab);
310  Range skin_edges;
311  CHKERR skin.find_skin(0, ents_1st_layer, false, skin_edges);
312  Range edges;
313  CHKERR moab.get_adjacencies(&*ents_1st_layer.begin(), 1, 1, false, edges);
314  EntityHandle meshset;
315  CHKERR mmanager_ptr->getMeshset(202, SIDESET, meshset);
316  CHKERR moab.add_entities(meshset, conn, 1);
317  CHKERR moab.add_entities(meshset, skin_edges);
318  }
319 
320  {
321  Range tets;
322  CHKERR moab.get_entities_by_type(0, MBTET, tets, false);
323  CHKERR m_field.partition_mesh(tets, 3, 2, n_partas);
324  }
325 
326  CHKERR moab.write_file("analysis_mesh.h5m");
327  }
328  CATCH_ERRORS;
329 
331 
332  return 0;
333 }
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
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]
constexpr int order
Calculate inverse of jacobian for face element.
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.
static char help[]
Definition: map_disp.cpp:33
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
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...
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.cpp.