52 PetscInt nb_sub_steps = 10;
53 PetscBool flg_file = PETSC_FALSE;
54 PetscBool is_partitioned = PETSC_FALSE;
55 PetscBool is_atom_test = PETSC_FALSE;
59 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
60 "Minimal surface area config",
"none");
61 CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
63 CHKERR PetscOptionsInt(
"-my_order",
"default approximation order",
"", 2,
65 CHKERR PetscOptionsInt(
"-my_nb_sub_steps",
"number of substeps",
"", 10,
66 &nb_sub_steps, PETSC_NULL);
67 CHKERR PetscOptionsBool(
"-my_is_partitioned",
68 "set if mesh is partitioned (this result that "
69 "each process keeps only part of the mesh",
70 "", PETSC_FALSE, &is_partitioned, PETSC_NULL);
73 "is used with testing, exit with error when diverged",
"",
74 PETSC_FALSE, &is_atom_test, PETSC_NULL);
76 ierr = PetscOptionsEnd();
80 if (flg_file != PETSC_TRUE) {
81 SETERRQ(PETSC_COMM_SELF, 1,
"*** ERROR -my_file (MESH FILE NEEDED)");
84 if (is_partitioned == PETSC_TRUE) {
87 option =
"PARALLEL=BCAST_DELETE;"
88 "PARALLEL_RESOLVE_SHARED_ENTS;"
89 "PARTITION=PARALLEL_PARTITION;";
96 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
98 pcomm =
new ParallelComm(&moab, PETSC_COMM_WORLD);
130 "MINIMAL_SURFACE_ELEMENT",
"U");
132 "MINIMAL_SURFACE_ELEMENT",
"U");
134 "MINIMAL_SURFACE_ELEMENT",
"U");
137 root_set, MBTRI,
"MINIMAL_SURFACE_ELEMENT");
144 Range bc_edges, bc_ents;
149 CHKERR moab.get_entities_by_type(root_set, MBTRI, tris,
false);
151 CHKERR skin.find_skin(root_set, tris,
false, bc_edges);
153 CHKERR moab.get_connectivity(bc_edges, bc_nodes);
154 bc_ents.merge(bc_edges);
155 bc_ents.merge(bc_nodes);
181 CHKERR DMCreate(PETSC_COMM_WORLD, &bc_dm);
182 CHKERR DMSetType(bc_dm,
"MoFEM");
187 CHKERR DMSetFromOptions(bc_dm);
194 CHKERR DMCreateGlobalVector(bc_dm, &T);
196 CHKERR DMCreateMatrix(bc_dm, &
A);
199 minimal_surface_element.feBcEdge.getOpPtrVector().push_back(
201 minimal_surface_element.feBcEdge.getOpPtrVector().push_back(
204 &minimal_surface_element.feBcEdge);
207 CHKERR MatAssemblyBegin(
A, MAT_FINAL_ASSEMBLY);
208 CHKERR MatAssemblyEnd(
A, MAT_FINAL_ASSEMBLY);
212 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
213 CHKERR KSPSetOperators(solver,
A,
A);
214 CHKERR KSPSetFromOptions(solver);
217 CHKERR KSPDestroy(&solver);
221 CHKERR VecGhostUpdateBegin(T, INSERT_VALUES, SCATTER_FORWARD);
222 CHKERR VecGhostUpdateEnd(T, INSERT_VALUES, SCATTER_FORWARD);
240 CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
241 CHKERR DMSetType(dm,
"MoFEM");
245 CHKERR DMSetFromOptions(dm);
256 CHKERR DMCreateGlobalVector(dm, &T);
266 auto det_ptr = boost::make_shared<VectorDouble>();
267 auto jac_ptr = boost::make_shared<MatrixDouble>();
268 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
269 minimal_surface_element.feRhs.getOpPtrVector().push_back(
271 minimal_surface_element.feRhs.getOpPtrVector().push_back(
274 minimal_surface_element.feRhs.getOpPtrVector().push_back(
277 minimal_surface_element.feRhs.getOpPtrVector().push_back(
279 minimal_surface_element.feRhs.getOpPtrVector().push_back(
281 "U", mse_common_data,
true));
282 minimal_surface_element.feRhs.getOpPtrVector().push_back(
287 minimal_surface_element.feLhs.getOpPtrVector().push_back(
289 minimal_surface_element.feLhs.getOpPtrVector().push_back(
292 minimal_surface_element.feLhs.getOpPtrVector().push_back(
295 minimal_surface_element.feLhs.getOpPtrVector().push_back(
297 minimal_surface_element.feLhs.getOpPtrVector().push_back(
299 "U", mse_common_data,
false));
300 minimal_surface_element.feLhs.getOpPtrVector().push_back(
307 &minimal_surface_element.feRhs, PETSC_NULL,
316 &minimal_surface_element.feLhs, PETSC_NULL,
322 SNESConvergedReason snes_reason;
327 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
332 CHKERR SNESSetFromOptions(snes);
335 CHKERR VecDuplicate(T, &T0);
338 double step_size = 1. / nb_sub_steps;
339 for (
int ss = 1; ss <= nb_sub_steps; ss++) {
340 CHKERR VecAXPY(T, step_size, T0);
342 CHKERR SNESSolve(snes, PETSC_NULL, T);
343 CHKERR SNESGetConvergedReason(snes, &snes_reason);
345 CHKERR SNESGetIterationNumber(snes, &its);
346 CHKERR PetscPrintf(PETSC_COMM_WORLD,
347 "number of Newton iterations = %D\n\n", its);
348 if (snes_reason < 0) {
351 "atom test diverged");
357 CHKERR VecGhostUpdateBegin(T, INSERT_VALUES, SCATTER_FORWARD);
358 CHKERR VecGhostUpdateEnd(T, INSERT_VALUES, SCATTER_FORWARD);
366 CHKERR SNESDestroy(&snes);
371 CHKERR post_proc.generateReferenceElementMesh();
372 CHKERR post_proc.addFieldValuesPostProc(
"U");
376 CHKERR post_proc.postProcMesh.write_file(
"out.h5m",
"MOAB",
377 "PARALLEL=WRITE_PART");