v0.15.0
Loading...
Searching...
No Matches
minimal_surface_area.cpp File Reference
#include <MoFEM.hpp>
#include <MethodForForceScaling.hpp>
#include <DirichletBC.hpp>
#include <PostProcOnRefMesh.hpp>
#include <MinimalSurfaceElement.hpp>

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 

Function Documentation

◆ main()

int main ( int argc,
char * argv[] )

Definition at line 41 of file minimal_surface_area.cpp.

41 {
42
43 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
44
45 try {
46
47 moab::Core mb_instance;
48 moab::Interface &moab = mb_instance;
49
50 char mesh_file_name[255];
51 PetscInt order = 2;
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",
62 mesh_file_name, 255, &flg_file);
63 CHKERR PetscOptionsInt("-my_order", "default approximation order", "", 2,
64 &order, PETSC_NULL);
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);
71 CHKERR PetscOptionsBool(
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 // Reade parameters from line command
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 // Read mesh to MOAB
85 const char *option;
86 option = "PARALLEL=BCAST_DELETE;"
87 "PARALLEL_RESOLVE_SHARED_ENTS;"
88 "PARTITION=PARALLEL_PARTITION;";
89 CHKERR moab.load_file(mesh_file_name, 0, option);
90 } else {
91 const char *option;
92 option = "";
93 CHKERR moab.load_file(mesh_file_name, 0, 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
100 MoFEM::Core core(moab);
101 MoFEM::Interface &m_field = core;
102 // meshset consisting all entities in mesh
103 EntityHandle root_set = moab.get_root_set();
104
105 // Seed all mesh entities to MoFEM database, those entities can be
106 // potentially used as finite elements or as entities which carry some
107 // approximation field.
108 BitRefLevel bit_level0;
109 bit_level0.set(0);
110 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
111 0, 2, bit_level0);
112
113 // Add field
114 CHKERR m_field.add_field("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
115 // Add entities to field (root_mesh, i.e. on all mesh entities fields are
116 // approx.) On those entities and lower dimension entities approximation is
117 // spaned,
118 CHKERR m_field.add_ents_to_field_by_type(root_set, MBTRI, "U");
119 // Set approximation oder
120 CHKERR m_field.set_field_order(root_set, MBTRI, "U", order);
121 CHKERR m_field.set_field_order(root_set, MBEDGE, "U", order);
122 CHKERR m_field.set_field_order(root_set, MBVERTEX, "U", 1);
123 CHKERR m_field.build_fields();
124
125 // Add finite element (this defines element declaration comes later)
126 CHKERR m_field.add_finite_element("MINIMAL_SURFACE_ELEMENT");
127 // Set on which fields that element operates
129 "MINIMAL_SURFACE_ELEMENT", "U");
131 "MINIMAL_SURFACE_ELEMENT", "U");
133 "MINIMAL_SURFACE_ELEMENT", "U");
134 // Add entities to that element
136 root_set, MBTRI, "MINIMAL_SURFACE_ELEMENT");
137
138 CHKERR m_field.add_finite_element("BC_ELEMENT");
139 CHKERR m_field.modify_finite_element_add_field_row("BC_ELEMENT", "U");
140 CHKERR m_field.modify_finite_element_add_field_col("BC_ELEMENT", "U");
141 CHKERR m_field.modify_finite_element_add_field_data("BC_ELEMENT", "U");
142 // Add entities to that element
143 Range bc_edges, bc_ents;
144 {
145 // Take a skin form the mesh
146 Skinner skin(&m_field.get_moab());
147 Range tris;
148 CHKERR moab.get_entities_by_type(root_set, MBTRI, tris, false);
149 Range bc_edges;
150 CHKERR skin.find_skin(root_set, tris, false, bc_edges);
151 Range bc_nodes;
152 CHKERR moab.get_connectivity(bc_edges, bc_nodes);
153 bc_ents.merge(bc_edges);
154 bc_ents.merge(bc_nodes);
155 // Add entities to that element
156 CHKERR m_field.add_ents_to_finite_element_by_type(bc_edges, MBEDGE,
157 "BC_ELEMENT");
158 }
159
160 // build finite elements
162 // build adjacencies between elements and degrees of freedom
163 CHKERR m_field.build_adjacencies(bit_level0);
164
165 // Create minimal surface area and elements
166 MinimalSurfaceElement minimal_surface_element(m_field);
167 // Data structure for common data
168 MinimalSurfaceElement::CommonData mse_common_data;
169
170 // This class is used to fix dofs on entities
171 DirichletFixFieldAtEntitiesBc fix_edges_ents(m_field, "U", NULL, NULL, NULL,
172 bc_ents);
173
174 CHKERR DMRegister_MoFEM("MoFEM");
175 // Solve BC problem
176 {
177 // define problems
178 DM bc_dm;
179 // Register DM problem
180 CHKERR DMCreate(PETSC_COMM_WORLD, &bc_dm);
181 CHKERR DMSetType(bc_dm, "MoFEM");
182 // Create DM instance
183 CHKERR DMMoFEMCreateMoFEM(bc_dm, &m_field, "DM_BC", bit_level0);
184 // Configure DM form line command options (DM itself, solvers,
185 // pre-conditioners, ... )
186 CHKERR DMSetFromOptions(bc_dm);
187 // Add elements to dm (only one here)
188 CHKERR DMMoFEMAddElement(bc_dm, "BC_ELEMENT");
189 // Set up problem (number dofs, partition mesh, etc.)
190 CHKERR DMSetUp(bc_dm);
191 Mat A;
192 Vec T, F;
193 CHKERR DMCreateGlobalVector(bc_dm, &T);
194 CHKERR VecDuplicate(T, &F);
195 CHKERR DMCreateMatrix(bc_dm, &A);
196
197 // Boundary problem
198 minimal_surface_element.feBcEdge.getOpPtrVector().push_back(
200 minimal_surface_element.feBcEdge.getOpPtrVector().push_back(
202 CHKERR DMoFEMLoopFiniteElements(bc_dm, "BC_ELEMENT",
203 &minimal_surface_element.feBcEdge);
204 CHKERR VecAssemblyBegin(F);
205 CHKERR VecAssemblyEnd(F);
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);
214 CHKERR KSPSetUp(solver);
215 CHKERR KSPSolve(solver, F, T);
216 CHKERR KSPDestroy(&solver);
217 }
218
219 // Scatter data on mesh
220 CHKERR VecGhostUpdateBegin(T, INSERT_VALUES, SCATTER_FORWARD);
221 CHKERR VecGhostUpdateEnd(T, INSERT_VALUES, SCATTER_FORWARD);
222 CHKERR DMoFEMMeshToGlobalVector(bc_dm, T, INSERT_VALUES, SCATTER_REVERSE);
223
224 CHKERR VecDestroy(&T);
225 CHKERR VecDestroy(&F);
226 CHKERR MatDestroy(&A);
227 CHKERR DMDestroy(&bc_dm);
228 }
229
230 // //define problems
231 // CHKERR m_field.add_problem("DM_MINIMAL_AREA");
232 // //set refinement level for problem
233 // CHKERR
234 // m_field.modify_problem_ref_level_add_bit("DM_MINIMAL_AREA",bit_level0);
235 // Create discrete manager instance
236 DM dm;
237 {
238 // Register problem
239 CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
240 CHKERR DMSetType(dm, "MoFEM");
241 // Create DM instance
242 CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "DM_MINIMAL_AREA", bit_level0);
243 // Get options fron line command
244 CHKERR DMSetFromOptions(dm);
245 // Add elements to dm (only one here)
246 CHKERR DMMoFEMAddElement(dm, "MINIMAL_SURFACE_ELEMENT");
247 // Set up data structures
248 CHKERR DMSetUp(dm);
249 }
250
251 // Create matrices and vectors used for analysis
252 Vec T, F;
253 Mat A;
254 {
255 CHKERR DMCreateGlobalVector(dm, &T);
256 CHKERR VecDuplicate(T, &F);
257 CHKERR DMCreateMatrix(dm, &A);
258 }
259
260 // Adding elements to DMSnes
261 {
262
263 // Right hand side operators (run for each entity in element)
264 // Calculate inverse of jacobian
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(
269 new OpCalculateHOJacForFace(jac_ptr));
270 minimal_surface_element.feRhs.getOpPtrVector().push_back(
271 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
272 // Apply inv_jac_ptr to shape functions
273 minimal_surface_element.feRhs.getOpPtrVector().push_back(
274 new OpSetInvJacH1ForFace(inv_jac_ptr));
275 // Below operators are specific to minimal area problem
276 minimal_surface_element.feRhs.getOpPtrVector().push_back(
277 new MinimalSurfaceElement::OpGetDataAtGaussPts("U", mse_common_data));
278 minimal_surface_element.feRhs.getOpPtrVector().push_back(
280 "U", mse_common_data, true));
281 minimal_surface_element.feRhs.getOpPtrVector().push_back(
282 new MinimalSurfaceElement::OpAssembleResidaul("U", mse_common_data));
283
284 // Left hand side operators (run for each entity in element)
285 // Calculate inverse of jacobian
286 minimal_surface_element.feLhs.getOpPtrVector().push_back(
287 new OpCalculateHOJacForFace(jac_ptr));
288 minimal_surface_element.feLhs.getOpPtrVector().push_back(
289 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
290 // Apply inv_jac_ptr to shape functions
291 minimal_surface_element.feLhs.getOpPtrVector().push_back(
292 new OpSetInvJacH1ForFace(inv_jac_ptr));
293 // Below operators are specific to minimal area problem
294 minimal_surface_element.feLhs.getOpPtrVector().push_back(
295 new MinimalSurfaceElement::OpGetDataAtGaussPts("U", mse_common_data));
296 minimal_surface_element.feLhs.getOpPtrVector().push_back(
298 "U", mse_common_data, false));
299 minimal_surface_element.feLhs.getOpPtrVector().push_back(
300 new MinimalSurfaceElement::OpAssembleTangent("U", mse_common_data));
301
302 // Rhs
303 CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, NULL, &fix_edges_ents,
304 NULL);
305 CHKERR DMMoFEMSNESSetFunction(dm, "MINIMAL_SURFACE_ELEMENT",
306 &minimal_surface_element.feRhs, PETSC_NULL,
307 PETSC_NULL);
309 &fix_edges_ents);
310
311 // Lhs
312 CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, NULL, &fix_edges_ents,
313 NULL);
314 CHKERR DMMoFEMSNESSetJacobian(dm, "MINIMAL_SURFACE_ELEMENT",
315 &minimal_surface_element.feLhs, PETSC_NULL,
316 PETSC_NULL);
318 &fix_edges_ents);
319 }
320
321 SNESConvergedReason snes_reason;
322 // Solve problem
323 {
324 SNES snes;
325 SnesCtx *snes_ctx;
326 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
327 // CHKERR SNESSetDM(snes,dm);
328 CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
329 CHKERR SNESSetFunction(snes, F, SnesRhs, snes_ctx);
330 CHKERR SNESSetJacobian(snes, A, A, SnesMat, snes_ctx);
331 CHKERR SNESSetFromOptions(snes);
332
333 Vec T0;
334 CHKERR VecDuplicate(T, &T0);
335 CHKERR DMoFEMMeshToLocalVector(dm, T0, INSERT_VALUES, SCATTER_FORWARD);
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);
340 CHKERR DMoFEMMeshToLocalVector(dm, T, INSERT_VALUES, SCATTER_REVERSE);
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) {
349 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
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);
358 CHKERR DMoFEMMeshToGlobalVector(dm, T, INSERT_VALUES, SCATTER_REVERSE);
359
360 // CHKERR MatView(A,PETSC_VIEWER_DRAW_SELF);
361 // std::string wait;
362 // std::cin >> wait;
363
364 CHKERR VecDestroy(&T0);
365 CHKERR SNESDestroy(&snes);
366 }
367
368 {
369 PostProcFaceOnRefinedMesh post_proc(m_field);
370 CHKERR post_proc.generateReferenceElementMesh();
371 CHKERR post_proc.addFieldValuesPostProc("U");
372 CHKERR DMoFEMLoopFiniteElements(dm, "MINIMAL_SURFACE_ELEMENT",
373 &post_proc);
374 // Save data on mesh
375 CHKERR post_proc.postProcMesh.write_file("out.h5m", "MOAB",
376 "PARALLEL=WRITE_PART");
377 }
378
379 // Clean and destroy
380 {
381 CHKERR DMDestroy(&dm);
382 CHKERR VecDestroy(&T);
383 CHKERR VecDestroy(&F);
384 CHKERR MatDestroy(&A);
385 }
386 }
388
390 return 0;
391}
#define DM_NO_ELEMENT
Definition DMMoFEM.hpp:10
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
@ H1
continuous field
Definition definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
#define CHKERR
Inline error check.
constexpr int order
@ F
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:488
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
Definition DMMoFEM.cpp:708
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 DMMoFEM.cpp:114
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition DMMoFEM.cpp:514
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
Definition DMMoFEM.cpp:749
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:576
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition DMMoFEM.cpp:1084
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition DMMoFEM.cpp:525
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.
char mesh_file_name[255]
static char help[]
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
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.
Definition SnesCtx.cpp:483
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.
Definition SnesCtx.cpp:223
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
constexpr AssemblyType A
Fix dofs on entities.
Implementation of minimal area element.
Managing BitRefLevels.
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.
Core (interface) class.
Definition Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:118
Deprecated interface functions.
Interface for nonlinear (SNES) solver.
Definition SnesCtx.hpp:15
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.

Variable Documentation

◆ help

char help[] = "...\n\n"
static

Definition at line 39 of file minimal_surface_area.cpp.