v0.9.0
minimal_surface_area.cpp
Go to the documentation of this file.
1 /** \file minimal_surface_area.cpp
2  \example minimal_surface_area.cpp
3  \brief Minimal surface area
4  \ingroup minimal_surface_area
5 
6  Implementation is based on:
7  <https://www.dealii.org/developer/doxygen/deal.II/step_15.html> Look for
8  formulation details on deal.II web page.
9 
10  \todo Make this example with heterogenous approx. apace
11  \todo Make this example with several refined meshes
12 
13 */
14 
15 /* This file is part of MoFEM.
16  * MoFEM is free software: you can redistribute it and/or modify it under
17  * the terms of the GNU Lesser General Public License as published by the
18  * Free Software Foundation, either version 3 of the License, or (at your
19  * option) any later version.
20  *
21  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
22  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
23  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
24  * License for more details.
25  *
26  * You should have received a copy of the GNU Lesser General Public
27  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
28 
29 #include <MoFEM.hpp>
30 using namespace MoFEM;
31 
33 #include <DirichletBC.hpp>
34 #include <PostProcOnRefMesh.hpp>
36 
37 using namespace MinimalSurfaceEquation;
38 
39 static char help[] = "...\n\n";
40 
41 int main(int argc, char *argv[]) {
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  CHKERR 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  ierr = PetscOptionsEnd();
77  CHKERRG(ierr);
78 
79  // Reade parameters from line command
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  // Read mesh to MOAB
86  const char *option;
87  option = "PARALLEL=BCAST_DELETE;"
88  "PARALLEL_RESOLVE_SHARED_ENTS;"
89  "PARTITION=PARALLEL_PARTITION;";
90  CHKERR moab.load_file(mesh_file_name, 0, option);
91  } else {
92  const char *option;
93  option = "";
94  CHKERR moab.load_file(mesh_file_name, 0, 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 
101  MoFEM::Core core(moab);
102  MoFEM::Interface &m_field = core;
103  // meshset consisting all entities in mesh
104  EntityHandle root_set = moab.get_root_set();
105 
106  // Seed all mesh entities to MoFEM database, those entities can be
107  // potentially used as finite elements or as entities which carry some
108  // approximation field.
109  BitRefLevel bit_level0;
110  bit_level0.set(0);
111  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
112  0, 2, bit_level0);
113 
114  // Add field
115  CHKERR m_field.add_field("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
116  // Add entities to field (root_mesh, i.e. on all mesh entities fields are
117  // approx.) On those entities and lower dimension entities approximation is
118  // spaned,
119  CHKERR m_field.add_ents_to_field_by_type(root_set, MBTRI, "U");
120  // Set approximation oder
121  CHKERR m_field.set_field_order(root_set, MBTRI, "U", order);
122  CHKERR m_field.set_field_order(root_set, MBEDGE, "U", order);
123  CHKERR m_field.set_field_order(root_set, MBVERTEX, "U", 1);
124  CHKERR m_field.build_fields();
125 
126  // Add finite element (this defines element declaration comes later)
127  CHKERR m_field.add_finite_element("MINIMAL_SURFACE_ELEMENT");
128  // Set on which fields that element operates
130  "MINIMAL_SURFACE_ELEMENT", "U");
132  "MINIMAL_SURFACE_ELEMENT", "U");
134  "MINIMAL_SURFACE_ELEMENT", "U");
135  // Add entities to that element
137  root_set, MBTRI, "MINIMAL_SURFACE_ELEMENT");
138 
139  CHKERR m_field.add_finite_element("BC_ELEMENT");
140  CHKERR m_field.modify_finite_element_add_field_row("BC_ELEMENT", "U");
141  CHKERR m_field.modify_finite_element_add_field_col("BC_ELEMENT", "U");
142  CHKERR m_field.modify_finite_element_add_field_data("BC_ELEMENT", "U");
143  // Add entities to that element
144  Range bc_edges, bc_ents;
145  {
146  // Take a skin form the mesh
147  Skinner skin(&m_field.get_moab());
148  Range tris;
149  CHKERR moab.get_entities_by_type(root_set, MBTRI, tris, false);
150  Range bc_edges;
151  CHKERR skin.find_skin(root_set, tris, false, bc_edges);
152  Range bc_nodes;
153  CHKERR moab.get_connectivity(bc_edges, bc_nodes);
154  bc_ents.merge(bc_edges);
155  bc_ents.merge(bc_nodes);
156  // Add entities to that element
157  CHKERR m_field.add_ents_to_finite_element_by_type(bc_edges, MBEDGE,
158  "BC_ELEMENT");
159  }
160 
161  // build finite elements
162  CHKERR m_field.build_finite_elements();
163  // build adjacencies between elements and degrees of freedom
164  CHKERR m_field.build_adjacencies(bit_level0);
165 
166  // Create minimal surface area and elements
167  MinimalSurfaceElement minimal_surface_element(m_field);
168  // Data structure for common data
169  MinimalSurfaceElement::CommonData mse_common_data;
170 
171  // This class is used to fix dofs on entities
172  DirichletFixFieldAtEntitiesBc fix_edges_ents(m_field, "U", NULL, NULL, NULL,
173  bc_ents);
174 
175  CHKERR DMRegister_MoFEM("MoFEM");
176  // Solve BC problem
177  {
178  // define problems
179  DM bc_dm;
180  // Register DM problem
181  CHKERR DMCreate(PETSC_COMM_WORLD, &bc_dm);
182  CHKERR DMSetType(bc_dm, "MoFEM");
183  // Create DM instance
184  CHKERR DMMoFEMCreateMoFEM(bc_dm, &m_field, "DM_BC", bit_level0);
185  // Configure DM form line command options (DM itself, solvers,
186  // pre-conditioners, ... )
187  CHKERR DMSetFromOptions(bc_dm);
188  // Add elements to dm (only one here)
189  CHKERR DMMoFEMAddElement(bc_dm, "BC_ELEMENT");
190  // Set up problem (number dofs, partition mesh, etc.)
191  CHKERR DMSetUp(bc_dm);
192  Mat A;
193  Vec T, F;
194  CHKERR DMCreateGlobalVector(bc_dm, &T);
195  CHKERR VecDuplicate(T, &F);
196  CHKERR DMCreateMatrix(bc_dm, &A);
197 
198  // Boundary problem
199  minimal_surface_element.feBcEdge.getOpPtrVector().push_back(
201  minimal_surface_element.feBcEdge.getOpPtrVector().push_back(
203  CHKERR DMoFEMLoopFiniteElements(bc_dm, "BC_ELEMENT",
204  &minimal_surface_element.feBcEdge);
205  CHKERR VecAssemblyBegin(F);
206  CHKERR VecAssemblyEnd(F);
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);
215  CHKERR KSPSetUp(solver);
216  CHKERR KSPSolve(solver, F, T);
217  CHKERR KSPDestroy(&solver);
218  }
219 
220  // Scatter data on mesh
221  CHKERR VecGhostUpdateBegin(T, INSERT_VALUES, SCATTER_FORWARD);
222  CHKERR VecGhostUpdateEnd(T, INSERT_VALUES, SCATTER_FORWARD);
223  CHKERR DMoFEMMeshToGlobalVector(bc_dm, T, INSERT_VALUES, SCATTER_REVERSE);
224 
225  CHKERR VecDestroy(&T);
226  CHKERR VecDestroy(&F);
227  CHKERR MatDestroy(&A);
228  CHKERR DMDestroy(&bc_dm);
229  }
230 
231  // //define problems
232  // CHKERR m_field.add_problem("DM_MINIMAL_AREA");
233  // //set refinement level for problem
234  // CHKERR
235  // m_field.modify_problem_ref_level_add_bit("DM_MINIMAL_AREA",bit_level0);
236  // Create discrete manager instance
237  DM dm;
238  {
239  // Register problem
240  CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
241  CHKERR DMSetType(dm, "MoFEM");
242  // Create DM instance
243  CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "DM_MINIMAL_AREA", bit_level0);
244  // Get options fron line command
245  CHKERR DMSetFromOptions(dm);
246  // Add elements to dm (only one here)
247  CHKERR DMMoFEMAddElement(dm, "MINIMAL_SURFACE_ELEMENT");
248  // Set up data structures
249  CHKERR DMSetUp(dm);
250  }
251 
252  // Create matrices and vectors used for analysis
253  Vec T, F;
254  Mat A;
255  {
256  CHKERR DMCreateGlobalVector(dm, &T);
257  CHKERR VecDuplicate(T, &F);
258  CHKERR DMCreateMatrix(dm, &A);
259  }
260 
261  // Adding elements to DMSnes
262  {
263 
264  // Right hand side operators (run for each entity in element)
265  // Calculate inverse of jacobian
266  minimal_surface_element.feRhs.getOpPtrVector().push_back(
267  new OpCalculateInvJacForFace(mse_common_data.invJac));
268  // Apply invJac to shape functions
269  minimal_surface_element.feRhs.getOpPtrVector().push_back(
270  new OpSetInvJacH1ForFace(mse_common_data.invJac));
271  // Below operators are specific to minimal area problem
272  minimal_surface_element.feRhs.getOpPtrVector().push_back(
273  new MinimalSurfaceElement::OpGetDataAtGaussPts("U", mse_common_data));
274  minimal_surface_element.feRhs.getOpPtrVector().push_back(
276  "U", mse_common_data, true));
277  minimal_surface_element.feRhs.getOpPtrVector().push_back(
278  new MinimalSurfaceElement::OpAssembleResidaul("U", mse_common_data));
279 
280  // Left hand side operators (run for each entity in element)
281  // Calculate inverse of jacobian
282  minimal_surface_element.feLhs.getOpPtrVector().push_back(
283  new OpCalculateInvJacForFace(mse_common_data.invJac));
284  // Apply invJac to shape functions
285  minimal_surface_element.feLhs.getOpPtrVector().push_back(
286  new OpSetInvJacH1ForFace(mse_common_data.invJac));
287  // Below operators are specific to minimal area problem
288  minimal_surface_element.feLhs.getOpPtrVector().push_back(
289  new MinimalSurfaceElement::OpGetDataAtGaussPts("U", mse_common_data));
290  minimal_surface_element.feLhs.getOpPtrVector().push_back(
292  "U", mse_common_data, false));
293  minimal_surface_element.feLhs.getOpPtrVector().push_back(
294  new MinimalSurfaceElement::OpAssembleTangent("U", mse_common_data));
295 
296  // Rhs
297  CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, NULL, &fix_edges_ents,
298  NULL);
299  CHKERR DMMoFEMSNESSetFunction(dm, "MINIMAL_SURFACE_ELEMENT",
300  &minimal_surface_element.feRhs, PETSC_NULL,
301  PETSC_NULL);
303  &fix_edges_ents);
304 
305  // Lhs
306  CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, NULL, &fix_edges_ents,
307  NULL);
308  CHKERR DMMoFEMSNESSetJacobian(dm, "MINIMAL_SURFACE_ELEMENT",
309  &minimal_surface_element.feLhs, PETSC_NULL,
310  PETSC_NULL);
312  &fix_edges_ents);
313  }
314 
315  SNESConvergedReason snes_reason;
316  // Solve problem
317  {
318  SNES snes;
319  SnesCtx *snes_ctx;
320  CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
321  // CHKERR SNESSetDM(snes,dm);
322  CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
323  CHKERR SNESSetFunction(snes, F, SnesRhs, snes_ctx);
324  CHKERR SNESSetJacobian(snes, A, A, SnesMat, snes_ctx);
325  CHKERR SNESSetFromOptions(snes);
326 
327  Vec T0;
328  CHKERR VecDuplicate(T, &T0);
329  CHKERR DMoFEMMeshToLocalVector(dm, T0, INSERT_VALUES, SCATTER_FORWARD);
330 
331  double step_size = 1. / nb_sub_steps;
332  for (int ss = 1; ss <= nb_sub_steps; ss++) {
333  CHKERR VecAXPY(T, step_size, T0);
334  CHKERR DMoFEMMeshToLocalVector(dm, T, INSERT_VALUES, SCATTER_REVERSE);
335  CHKERR SNESSolve(snes, PETSC_NULL, T);
336  CHKERR SNESGetConvergedReason(snes, &snes_reason);
337  int its;
338  CHKERR SNESGetIterationNumber(snes, &its);
339  CHKERR PetscPrintf(PETSC_COMM_WORLD,
340  "number of Newton iterations = %D\n\n", its);
341  if (snes_reason < 0) {
342  if (is_atom_test) {
343  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
344  "atom test diverged");
345  } else {
346  break;
347  }
348  }
349  }
350  CHKERR VecGhostUpdateBegin(T, INSERT_VALUES, SCATTER_FORWARD);
351  CHKERR VecGhostUpdateEnd(T, INSERT_VALUES, SCATTER_FORWARD);
352  CHKERR DMoFEMMeshToGlobalVector(dm, T, INSERT_VALUES, SCATTER_REVERSE);
353 
354  // CHKERR MatView(A,PETSC_VIEWER_DRAW_SELF);
355  // std::string wait;
356  // std::cin >> wait;
357 
358  CHKERR VecDestroy(&T0);
359  CHKERR SNESDestroy(&snes);
360  }
361 
362  {
363  PostProcFaceOnRefinedMesh post_proc(m_field);
365  CHKERR post_proc.addFieldValuesPostProc("U");
366  CHKERR DMoFEMLoopFiniteElements(dm, "MINIMAL_SURFACE_ELEMENT",
367  &post_proc);
368  // Save data on mesh
369  CHKERR post_proc.postProcMesh.write_file("out.h5m", "MOAB",
370  "PARALLEL=WRITE_PART");
371  }
372 
373  // Clean and destroy
374  {
375  CHKERR DMDestroy(&dm);
376  CHKERR VecDestroy(&T);
377  CHKERR VecDestroy(&F);
378  CHKERR MatDestroy(&A);
379  }
380  }
381  CATCH_ERRORS;
382 
384  return 0;
385 }
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
EdgeElement feBcEdge
Used to calculate dofs on boundary.
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:22
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: DMMMoFEM.cpp:626
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 build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
Implementation of minimal surface elementImplementation is based on: https://www.dealii....
int main(int argc, char *argv[])
virtual moab::Interface & get_moab()=0
Post-process fields on refined mesh.
moab::Interface & postProcMesh
Implementation of minimal area element.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:434
SurfaceElement feLhs
To calculate left hand side.
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:414
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:121
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:544
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.
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.
Core (interface) class.
Definition: Core.hpp:50
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:496
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
char mesh_file_name[255]
static char help[]
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:51
MoFEMErrorCode generateReferenceElementMesh()
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
Calculate inverse of jacobian for face element.
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMMoFEM.cpp:942
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.
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:446
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:27
Managing BitRefLevels.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:53
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:146
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: DMMMoFEM.cpp:667
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: DMMMoFEM.cpp:109
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
#define CHKERR
Inline error check.
Definition: definitions.h:596
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
Transform local reference derivatives of shape functions to global derivatives.
Fix dofs on entities.
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:285
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
Postprocess on face.
Evaluate function values and gradients at Gauss Pts.
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
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:433
continuous field
Definition: definitions.h:171
SurfaceElement feRhs
To calculate right hand side.
Minimal surface equation.
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:23
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:61