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 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);
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 ierr = PetscOptionsEnd();
78
79
80 if (flg_file != PETSC_TRUE) {
81 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
82 }
83
84 if (is_partitioned == PETSC_TRUE) {
85
86 const char *option;
87 option = "PARALLEL=BCAST_DELETE;"
88 "PARALLEL_RESOLVE_SHARED_ENTS;"
89 "PARTITION=PARALLEL_PARTITION;";
91 } else {
92 const char *option;
93 option = "";
95 }
96 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
97 if (pcomm == NULL)
98 pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
99 }
100
103
105
106
107
108
110 bit_level0.set(0);
112 0, 2, bit_level0);
113
114
116
117
118
120
125
126
128
130 "MINIMAL_SURFACE_ELEMENT", "U");
132 "MINIMAL_SURFACE_ELEMENT", "U");
134 "MINIMAL_SURFACE_ELEMENT", "U");
135
137 root_set, MBTRI, "MINIMAL_SURFACE_ELEMENT");
138
143
144 Range bc_edges, bc_ents;
145 {
146
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);
156
158 "BC_ELEMENT");
159 }
160
161
163
165
166
168
170
171
173 bc_ents);
174
176
177 {
178
179 DM bc_dm;
180
181 CHKERR DMCreate(PETSC_COMM_WORLD, &bc_dm);
182 CHKERR DMSetType(bc_dm,
"MoFEM");
183
185
186
187 CHKERR DMSetFromOptions(bc_dm);
188
190
194 CHKERR DMCreateGlobalVector(bc_dm, &T);
196 CHKERR DMCreateMatrix(bc_dm, &A);
197
198
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);
209
210 {
211 KSP solver;
212 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
213 CHKERR KSPSetOperators(solver, A, A);
214 CHKERR KSPSetFromOptions(solver);
217 CHKERR KSPDestroy(&solver);
218 }
219
220
221 CHKERR VecGhostUpdateBegin(T, INSERT_VALUES, SCATTER_FORWARD);
222 CHKERR VecGhostUpdateEnd(T, INSERT_VALUES, SCATTER_FORWARD);
224
229 }
230
231
232
233
234
235
236
237 DM dm;
238 {
239
240 CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
241 CHKERR DMSetType(dm,
"MoFEM");
242
244
245 CHKERR DMSetFromOptions(dm);
246
248
250 }
251
252
255 {
256 CHKERR DMCreateGlobalVector(dm, &T);
258 CHKERR DMCreateMatrix(dm, &A);
259 }
260
261
262 {
263
264
265
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(
273
274 minimal_surface_element.feRhs.getOpPtrVector().push_back(
276
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(
284
285
286
287 minimal_surface_element.feLhs.getOpPtrVector().push_back(
289 minimal_surface_element.feLhs.getOpPtrVector().push_back(
291
292 minimal_surface_element.feLhs.getOpPtrVector().push_back(
294
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(
302
303
305 NULL);
307 &minimal_surface_element.feRhs, PETSC_NULL,
308 PETSC_NULL);
310 &fix_edges_ents);
311
312
314 NULL);
316 &minimal_surface_element.feLhs, PETSC_NULL,
317 PETSC_NULL);
319 &fix_edges_ents);
320 }
321
322 SNESConvergedReason snes_reason;
323
324 {
325 SNES snes;
327 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
328
332 CHKERR SNESSetFromOptions(snes);
333
335 CHKERR VecDuplicate(T, &T0);
337
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);
344 int its;
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) {
349 if (is_atom_test) {
351 "atom test diverged");
352 } else {
353 break;
354 }
355 }
356 }
357 CHKERR VecGhostUpdateBegin(T, INSERT_VALUES, SCATTER_FORWARD);
358 CHKERR VecGhostUpdateEnd(T, INSERT_VALUES, SCATTER_FORWARD);
360
361
362
363
364
366 CHKERR SNESDestroy(&snes);
367 }
368
369 {
371 CHKERR post_proc.generateReferenceElementMesh();
372 CHKERR post_proc.addFieldValuesPostProc(
"U");
374 &post_proc);
375
376 CHKERR post_proc.postProcMesh.write_file(
"out.h5m",
"MOAB",
377 "PARALLEL=WRITE_PART");
378 }
379
380
381 {
386 }
387 }
389
391 return 0;
392}
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MYPCOMM_INDEX
default communicator number PCOMM
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ 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_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
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 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.
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
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 refernce to pointer of interface.