18 using namespace MoFEM;
26 #include <boost/tokenizer.hpp>
30 using namespace boost;
35 int main(
int argc,
char *argv[]) {
45 char data_file_name1[255];
46 char data_file_name2[255];
47 PetscBool flg_file = PETSC_FALSE;
49 PetscBool flg_n_part = PETSC_FALSE;
50 PetscInt n_partas = 1;
51 PetscBool flg_thickness = PETSC_FALSE;
54 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"none",
"none");
55 CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
58 CHKERR PetscOptionsString(
"-my_data_x",
"data file name",
"",
"data_x.h5m",
59 data_file_name1, 255, PETSC_NULL);
61 CHKERR PetscOptionsString(
"-my_data_y",
"data file name",
"",
"data_y.h5m",
62 data_file_name2, 255, PETSC_NULL);
64 CHKERR PetscOptionsInt(
"-my_order",
"default approximation order",
"", 1,
67 CHKERR PetscOptionsInt(
"-my_nparts",
"number of parts",
"", 1, &n_partas,
70 CHKERR PetscOptionsScalar(
"-my_thinckness",
"top layer thickness",
"",
71 thickness, &thickness, &flg_thickness);
73 ierr = PetscOptionsEnd();
76 if (flg_file != PETSC_TRUE) {
77 SETERRQ(PETSC_COMM_SELF, 1,
"*** ERROR -my_file (MESH FILE NEEDED)");
80 if (flg_n_part != PETSC_TRUE) {
81 SETERRQ(PETSC_COMM_SELF, 1,
"*** ERROR partitioning number not given");
84 if (flg_thickness != PETSC_TRUE) {
85 SETERRQ(PETSC_COMM_SELF, 1,
"*** ERROR thickness of top layer not given");
91 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
93 pcomm =
new ParallelComm(&moab, PETSC_COMM_WORLD);
150 DMType dm_name =
"DM_DISP";
153 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_disp);
154 CHKERR DMSetType(dm_disp, dm_name);
159 CHKERR DMSetFromOptions(dm_disp);
170 string file1(data_file_name1);
175 string file2(data_file_name2);
180 boost::ptr_vector<MoFEM::ForcesAndSourcesCore::UserDataOperator>
185 CHKERR DMCreateMatrix(dm_disp, &
A);
191 CHKERR VecDuplicate(Fx, &Fy);
192 CHKERR VecDuplicate(Dx, &Dy);
195 CHKERR VecZeroEntries(Fx);
196 CHKERR VecZeroEntries(Fy);
201 common_data.
globF = Fx;
202 list_of_operators.push_back(
203 new OpCalculateInvJacForFace(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(
209 list_of_operators.push_back(
new OpAssmbleRhs(
"DISP_X", common_data));
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);
222 common_data.
globF = Fy;
223 list_of_operators.clear();
224 list_of_operators.push_back(
225 new OpCalculateInvJacForFace(common_data.
invJac));
227 list_of_operators.push_back(
229 list_of_operators.push_back(
new OpAssmbleRhs(
"DISP_X", common_data));
234 CHKERR VecGhostUpdateBegin(Fy, ADD_VALUES, SCATTER_REVERSE);
235 CHKERR VecGhostUpdateEnd(Fy, ADD_VALUES, SCATTER_REVERSE);
236 CHKERR VecAssemblyBegin(Fy);
237 CHKERR VecAssemblyEnd(Fy);
247 CHKERR KSPCreate(PETSC_COMM_WORLD, &ksp);
249 CHKERR KSPSetFromOptions(ksp);
250 CHKERR KSPSolve(ksp, Fx, Dx);
251 CHKERR KSPSolve(ksp, Fy, Dy);
254 CHKERR VecGhostUpdateBegin(Dx, INSERT_VALUES, SCATTER_FORWARD);
255 CHKERR VecGhostUpdateEnd(Dx, INSERT_VALUES, SCATTER_FORWARD);
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);
261 "DM_DISP",
"DISP_X",
"DISP_Y",
COL, Dy, INSERT_VALUES, SCATTER_REVERSE);
268 sum_dx += (*dof)->getFieldData();
274 (*dof)->getFieldData() -= sum_dx;
280 sum_dy += (*dof)->getFieldData();
286 (*dof)->getFieldData() -= sum_dy;
296 "PARALLEL=WRITE_PART");
303 CHKERR DMDestroy(&dm_disp);
312 Range ents_2nd_layer;
314 ents_2nd_layer,
true);
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();
321 CHKERR moab.get_coords(&*nit, 1, coords);
322 coords[2] += thickness;
323 CHKERR moab.create_vertex(coords, map_nodes[*nit]);
328 for (Range::iterator tit = ents_2nd_layer.begin();
329 tit != ents_2nd_layer.end(); tit++) {
332 CHKERR moab.get_connectivity(*tit, conn, num_nodes);
334 for (
int n = 0;
n != 3;
n++) {
335 conn_prism[
n] = conn[
n];
336 conn_prism[
n + 3] = map_nodes[conn[
n]];
339 CHKERR moab.create_element(MBPRISM, conn_prism, 6, prism);
341 CHKERR moab.get_adjacencies(&prism, 1, 2,
true, adj);
342 CHKERR moab.get_adjacencies(&prism, 1, 1,
true, adj);
343 prisms.insert(prism);
345 CHKERR moab.side_element(prism, 2, 4, tri);
346 top_tris.insert(tri);
353 mat_elastic.
data.Young = 1;
354 mat_elastic.
data.Poisson = 0;
363 Range ents_1st_layer;
365 ents_1st_layer,
true);
368 CHKERR moab.get_connectivity(ents_1st_layer[0], conn, num_nodes);
371 CHKERR skin.find_skin(0, ents_1st_layer,
false, skin_edges);
373 CHKERR moab.get_adjacencies(&*ents_1st_layer.begin(), 1, 1,
false, edges);
376 CHKERR moab.add_entities(meshset, conn, 1);
377 CHKERR moab.add_entities(meshset, skin_edges);
382 CHKERR moab.get_entities_by_dimension(0, 3, ents3d,
false);
386 CHKERR moab.write_file(
"analysis_mesh.h5m");
392 PetscFunctionReturn(0);