41 {
42
44
45 try {
46
47 moab::Core mb_instance;
48 moab::Interface &moab = mb_instance;
49
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;
56
57 {
58
59 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);
72 "-my_is_atom_test",
73 "is used with testing, exit with error when diverged", "",
74 PETSC_FALSE, &is_atom_test, PETSC_NULL);
75
76 PetscOptionsEnd();
77
78
79 if (flg_file != PETSC_TRUE) {
80 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
81 }
82
83 if (is_partitioned == PETSC_TRUE) {
84
85 const char *option;
86 option = "PARALLEL=BCAST_DELETE;"
87 "PARALLEL_RESOLVE_SHARED_ENTS;"
88 "PARTITION=PARALLEL_PARTITION;";
90 } else {
91 const char *option;
92 option = "";
94 }
95 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
96 if (pcomm == NULL)
97 pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
98 }
99
102
104
105
106
107
109 bit_level0.set(0);
111 0, 2, bit_level0);
112
113
115
116
117
119
124
125
127
129 "MINIMAL_SURFACE_ELEMENT", "U");
131 "MINIMAL_SURFACE_ELEMENT", "U");
133 "MINIMAL_SURFACE_ELEMENT", "U");
134
136 root_set, MBTRI, "MINIMAL_SURFACE_ELEMENT");
137
142
143 Range bc_edges, bc_ents;
144 {
145
148 CHKERR moab.get_entities_by_type(root_set, MBTRI, tris,
false);
150 CHKERR skin.find_skin(root_set, tris,
false, bc_edges);
152 CHKERR moab.get_connectivity(bc_edges, bc_nodes);
153 bc_ents.merge(bc_edges);
154 bc_ents.merge(bc_nodes);
155
157 "BC_ELEMENT");
158 }
159
160
162
164
165
167
169
170
172 bc_ents);
173
175
176 {
177
178 DM bc_dm;
179
180 CHKERR DMCreate(PETSC_COMM_WORLD, &bc_dm);
181 CHKERR DMSetType(bc_dm,
"MoFEM");
182
184
185
186 CHKERR DMSetFromOptions(bc_dm);
187
189
193 CHKERR DMCreateGlobalVector(bc_dm, &T);
195 CHKERR DMCreateMatrix(bc_dm, &A);
196
197
198 minimal_surface_element.feBcEdge.getOpPtrVector().push_back(
200 minimal_surface_element.feBcEdge.getOpPtrVector().push_back(
203 &minimal_surface_element.feBcEdge);
206 CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
207 CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
208
209 {
210 KSP solver;
211 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
212 CHKERR KSPSetOperators(solver, A, A);
213 CHKERR KSPSetFromOptions(solver);
216 CHKERR KSPDestroy(&solver);
217 }
218
219
220 CHKERR VecGhostUpdateBegin(T, INSERT_VALUES, SCATTER_FORWARD);
221 CHKERR VecGhostUpdateEnd(T, INSERT_VALUES, SCATTER_FORWARD);
223
228 }
229
230
231
232
233
234
235
236 DM dm;
237 {
238
239 CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
240 CHKERR DMSetType(dm,
"MoFEM");
241
243
244 CHKERR DMSetFromOptions(dm);
245
247
249 }
250
251
254 {
255 CHKERR DMCreateGlobalVector(dm, &T);
257 CHKERR DMCreateMatrix(dm, &A);
258 }
259
260
261 {
262
263
264
265 auto det_ptr = boost::make_shared<VectorDouble>();
266 auto jac_ptr = boost::make_shared<MatrixDouble>();
267 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
268 minimal_surface_element.feRhs.getOpPtrVector().push_back(
270 minimal_surface_element.feRhs.getOpPtrVector().push_back(
272
273 minimal_surface_element.feRhs.getOpPtrVector().push_back(
275
276 minimal_surface_element.feRhs.getOpPtrVector().push_back(
278 minimal_surface_element.feRhs.getOpPtrVector().push_back(
280 "U", mse_common_data, true));
281 minimal_surface_element.feRhs.getOpPtrVector().push_back(
283
284
285
286 minimal_surface_element.feLhs.getOpPtrVector().push_back(
288 minimal_surface_element.feLhs.getOpPtrVector().push_back(
290
291 minimal_surface_element.feLhs.getOpPtrVector().push_back(
293
294 minimal_surface_element.feLhs.getOpPtrVector().push_back(
296 minimal_surface_element.feLhs.getOpPtrVector().push_back(
298 "U", mse_common_data, false));
299 minimal_surface_element.feLhs.getOpPtrVector().push_back(
301
302
304 NULL);
306 &minimal_surface_element.feRhs, PETSC_NULL,
307 PETSC_NULL);
309 &fix_edges_ents);
310
311
313 NULL);
315 &minimal_surface_element.feLhs, PETSC_NULL,
316 PETSC_NULL);
318 &fix_edges_ents);
319 }
320
321 SNESConvergedReason snes_reason;
322
323 {
324 SNES snes;
326 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
327
331 CHKERR SNESSetFromOptions(snes);
332
333 Vec T0;
334 CHKERR VecDuplicate(T, &T0);
336
337 double step_size = 1. / nb_sub_steps;
338 for (int ss = 1; ss <= nb_sub_steps; ss++) {
339 CHKERR VecAXPY(T, step_size, T0);
341 CHKERR SNESSolve(snes, PETSC_NULL, T);
342 CHKERR SNESGetConvergedReason(snes, &snes_reason);
343 int its;
344 CHKERR SNESGetIterationNumber(snes, &its);
345 CHKERR PetscPrintf(PETSC_COMM_WORLD,
346 "number of Newton iterations = %D\n\n", its);
347 if (snes_reason < 0) {
348 if (is_atom_test) {
350 "atom test diverged");
351 } else {
352 break;
353 }
354 }
355 }
356 CHKERR VecGhostUpdateBegin(T, INSERT_VALUES, SCATTER_FORWARD);
357 CHKERR VecGhostUpdateEnd(T, INSERT_VALUES, SCATTER_FORWARD);
359
360
361
362
363
365 CHKERR SNESDestroy(&snes);
366 }
367
368 {
370 CHKERR post_proc.generateReferenceElementMesh();
371 CHKERR post_proc.addFieldValuesPostProc(
"U");
373 &post_proc);
374
375 CHKERR post_proc.postProcMesh.write_file(
"out.h5m",
"MOAB",
376 "PARALLEL=WRITE_PART");
377 }
378
379
380 {
385 }
386 }
388
390 return 0;
391}
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
#define MYPCOMM_INDEX
default communicator number PCOMM
@ MOFEM_ATOM_TEST_INVALID
#define CHKERR
Inline error check.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
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.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
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_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
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.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
Keep date between operators.
Integrate vector on lhs,.
Integrate vector on rhs,.
Evaluate function values and gradients at Gauss Pts.
Evaluate function values and gradients at Gauss Pts.
Implementation of minimal area element.
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Interface for nonlinear (SNES) solver.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.