v0.9.0
MinimalSurfaceElement.hpp
Go to the documentation of this file.
1 /** \file MinimalSurfaceElement.hpp
2  \brief Implementation of minimal surface element
3  \ingroup minimal_surface_area
4 
5  Implementation is based on:
6  <https://www.dealii.org/developer/doxygen/deal.II/step_15.html>
7 
8 */
9 
10 /* This file is part of MoFEM.
11  * MoFEM is free software: you can redistribute it and/or modify it under
12  * the terms of the GNU Lesser General Public License as published by the
13  * Free Software Foundation, either version 3 of the License, or (at your
14  * option) any later version.
15  *
16  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
17  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19  * License for more details.
20  *
21  * You should have received a copy of the GNU Lesser General Public
22  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
23 
24 #ifndef __MINIMAL_SURFACE_ELEMENT_HPP__
25 #define __MINIMAL_SURFACE_ELEMENT_HPP__
26 
27 /** \brief Minimal surface equation
28 \ingroup minimal_surface_area
29 
30 */
32 
33 /** \brief Implementation of minimal area element
34 \ingroup minimal_surface_area
35 
36 */
38 
40 
44 
45  /** \brief Set integrate rule
46  */
47  int getRule(int order) { return 2 * (order - 1) + 1; };
48  };
49 
50  /** \brief Bondary problem
51 
52  \f[
53  \int_L \mathbf{N}^\textrm{T}\mathbf{N} \textrm{d}L \; \mathbf{u} = \int_L
54  \mathbf{N}^\textrm{T} g(x,y) \textrm{d}L \f]
55 
56  */
60 
61  /** \brief Set integrate rule
62  */
63  int getRule(int order) { return 2 * order; };
64  };
65 
66  SurfaceElement feRhs; ///< To calculate right hand side
67  SurfaceElement feLhs; ///< To calculate left hand side
68  EdgeElement feBcEdge; ///< Used to calculate dofs on boundary
69 
71  : mField(m_field), feRhs(m_field), feLhs(m_field), feBcEdge(m_field) {}
72 
73  /** \brief Keep date between operators
74  */
75  struct CommonData {
76  MatrixDouble gradU; ///< nb_gauss_pts x 2 gradients at integration pts
77  VectorDouble normGradU2; ///< size of nb_gauss_pts, norm of gradient
79  aN; ///< size of nb_gauss_pts, \f$a_n = \frac{1}{\sqrt{1+\|u_i\|^2}}\f$
80  MatrixDouble aNbyGradU; ///< nb_gauss_pts x 2, \f$a_n \nabla u\f$
81  MatrixDouble aNpow3byGradU; ///< nb_gauss_pts x 2, \f$a_n^3 \nabla u\f$
82  MatrixDouble invJac; ///< inverse of jacobian matrix
83  };
84 
85  /** \brief Integrate vector on rhs,
86 
87  \f[
88  \mathbf{F} = \int_L \mathbf{N}^\textrm{T} g(x,y) \textrm{d}L
89  \f]
90 
91  */
94  Vec vF;
95  OpAssmebleBcRhs(const string field_name, Vec v_f)
97  field_name, UserDataOperator::OPROW),
98  vF(v_f) {}
99 
100  /** \brief Function on boundary
101  Inherit this class and overload this function to change bc.
102  */
103  virtual double evalFunction(double x, double y) {
104  return sin(2 * M_PI * (x + y));
105  // return sin(M_PI*x);
106  }
107 
109  PetscErrorCode doWork(int side, EntityType type,
112 
113  int nb_dofs = data.getFieldData().size();
114  if (nb_dofs == 0) {
116  }
117  nF.resize(nb_dofs, false);
118  nF.clear();
119  int nb_gauss_pts = data.getN().size1();
120  for (int gg = 0; gg != nb_gauss_pts; gg++) {
121  double val = getLength() * getGaussPts()(1, gg);
122  double x = getCoordsAtGaussPts()(gg, 0);
123  double y = getCoordsAtGaussPts()(gg, 1);
124  noalias(nF) += val * evalFunction(x, y) * data.getN(gg);
125  }
126  CHKERR VecSetValues(vF, data.getIndices().size(),
127  &*data.getIndices().data().begin(), &nF[0],
128  ADD_VALUES);
130  }
131  };
132 
133  /** \brief Integrate vector on lhs,
134 
135  \f[
136  \mathbf{A} = \int_L \mathbf{N}^\textrm{T}\mathbf{N} \textrm{d}L
137  \f]
138 
139  */
142  Mat mA;
143  OpAssmebleBcLhs(const string field_name, Mat m_a)
145  field_name, UserDataOperator::OPROWCOL),
146  mA(m_a) {}
148  PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
149  EntityType col_type,
153 
154  int row_nb_dofs = row_data.getIndices().size();
155  if (row_nb_dofs == 0) {
157  }
158  int col_nb_dofs = col_data.getIndices().size();
159  if (col_nb_dofs == 0) {
161  }
162  K.resize(row_nb_dofs, col_nb_dofs, false);
163  K.clear();
164  int nb_gauss_pts = row_data.getN().size1();
165  for (int gg = 0; gg != nb_gauss_pts; gg++) {
166  double val = getLength() * getGaussPts()(1, gg);
167  VectorAdaptor row_n = row_data.getN(gg);
168  VectorAdaptor col_n = col_data.getN(gg);
169  noalias(K) += val * outer_prod(row_n, col_n);
170  }
171  CHKERR MatSetValues(mA, row_nb_dofs,
172  &*row_data.getIndices().data().begin(), col_nb_dofs,
173  &*col_data.getIndices().data().begin(),
174  &*K.data().begin(), ADD_VALUES);
175  // Matrix is symmetric, assemble off diagonal part
176  if (row_type != col_type || row_side != col_side) {
177  transK.resize(col_nb_dofs, row_nb_dofs, false);
178  noalias(transK) = trans(K);
179  CHKERR MatSetValues(mA, col_nb_dofs,
180  &*col_data.getIndices().data().begin(), row_nb_dofs,
181  &*row_data.getIndices().data().begin(),
182  &*transK.data().begin(), ADD_VALUES);
183  }
184 
186  }
187  };
188 
189  /** \brief Evaluate function values and gradients at Gauss Pts
190  */
194  OpGetDataAtGaussPts(const string field_name, CommonData &common_data)
196  field_name, UserDataOperator::OPCOL),
197  commonData(common_data) {}
198 
199  PetscErrorCode doWork(int side, EntityType type,
202 
203  int nb_dofs = data.getFieldData().size();
204  if (nb_dofs == 0) {
206  }
207  int nb_gauss_pts = data.getN().size1();
208 
209  // Element loops over entities, it start from vertices
210  if (type == MBVERTEX) {
211  // Setting size of common data vectors
212  commonData.gradU.resize(nb_gauss_pts, 2);
213  commonData.gradU.clear();
214  }
215 
216  for (int gg = 0; gg != nb_gauss_pts; gg++) {
217  MatrixAdaptor grad_shape_fun = data.getDiffN(gg);
218  ublas::matrix_row<MatrixDouble> grad_at_gauss_pt(commonData.gradU, gg);
219  noalias(grad_at_gauss_pt) +=
220  prod(trans(grad_shape_fun), data.getFieldData());
221  }
222 
224  }
225  };
226 
227  /** \brief Evaluate function values and gradients at Gauss Pts
228  */
233  OpCalculateCoefficientsAtGaussPts(const string field_name,
234  CommonData &common_data,
235  bool rhs_and_not_lhs)
237  field_name, UserDataOperator::OPCOL),
238  commonData(common_data), rhsAndNotLhs(rhs_and_not_lhs) {}
239 
240  PetscErrorCode doWork(int side, EntityType type,
243 
244  // Element loops over entities, it start from vertices
245  if (type == MBVERTEX) {
246  int nb_gauss_pts = data.getN().size1();
247  commonData.normGradU2.resize(nb_gauss_pts, false);
248  for (int gg = 0; gg != nb_gauss_pts; gg++) {
249  ublas::matrix_row<MatrixDouble> grad_at_gauss_pt(commonData.gradU,
250  gg);
251  commonData.normGradU2[gg] =
252  inner_prod(grad_at_gauss_pt, grad_at_gauss_pt);
253  }
254  commonData.aN.resize(nb_gauss_pts, false);
255  for (int gg = 0; gg != nb_gauss_pts; gg++) {
256  commonData.aN[gg] = 1. / sqrt(1 + commonData.normGradU2[gg]);
257  }
258  if (rhsAndNotLhs) {
259  // needed at right hand side when residual is calculated
260  commonData.aNbyGradU.resize(nb_gauss_pts, 2, false);
261  for (int gg = 0; gg != nb_gauss_pts; gg++) {
262  ublas::matrix_row<MatrixDouble> grad_at_gauss_pt(commonData.gradU,
263  gg);
264  for (int rr = 0; rr != 2; rr++) {
265  commonData.aNbyGradU(gg, rr) =
266  commonData.aN[gg] * grad_at_gauss_pt[rr];
267  }
268  }
269  } else {
270  // needed at left hand side when matrix is calculated
271  commonData.aNpow3byGradU.resize(nb_gauss_pts, 2, false);
272  for (int gg = 0; gg != nb_gauss_pts; gg++) {
273  ublas::matrix_row<MatrixDouble> grad_at_gauss_pt(commonData.gradU,
274  gg);
275  for (int rr = 0; rr != 2; rr++) {
276  commonData.aNpow3byGradU(gg, rr) =
277  pow(commonData.aN[gg], 3) * grad_at_gauss_pt[rr];
278  }
279  }
280  }
281  }
282 
283  // doVerticesRow = true;
284  // doEdgesRow = false;
285  // doQuadsRow = false;
286  // doTrisRow = false;
287  // doTetsRow = false;
288  // doPrismsRow = false;
289 
291  }
292  };
293 
294  /** \brief Assemble residual
295  */
299  OpAssembleResidaul(const string field_name, CommonData &common_data)
301  field_name, UserDataOperator::OPROW),
302  commonData(common_data) {}
304  PetscErrorCode doWork(int side, EntityType type,
307 
308  int nb_dofs = data.getIndices().size();
309  if (nb_dofs == 0) {
311  }
312  rEsidual.resize(nb_dofs, false);
313  rEsidual.clear();
314  int nb_gauss_pts = data.getN().size1();
315  for (int gg = 0; gg != nb_gauss_pts; gg++) {
316  double val = getGaussPts()(2, gg) * getArea();
317  ublas::matrix_row<MatrixDouble> an_by_grad_u(commonData.aNbyGradU, gg);
318  noalias(rEsidual) += val * prod(data.getDiffN(gg), an_by_grad_u);
319  }
320  CHKERR VecSetValues(getFEMethod()->snes_f, data.getIndices().size(),
321  &*data.getIndices().data().begin(), &rEsidual[0],
322  ADD_VALUES);
324  }
325  };
326 
327  /** \brief Assemble tangent matrix
328  */
332  OpAssembleTangent(const string field_name, CommonData &common_data)
334  field_name, UserDataOperator::OPROWCOL),
335  commonData(common_data) {
336  sYmm = false; // matrix is not symmetric
337  }
341  PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
342  EntityType col_type,
346 
347  int row_nb_dofs = row_data.getIndices().size();
348  if (row_nb_dofs == 0) {
350  }
351  int col_nb_dofs = col_data.getIndices().size();
352  if (col_nb_dofs == 0) {
354  }
355  kTangent.resize(row_nb_dofs, col_nb_dofs, false);
356  kTangent.clear();
357  auxVec.resize(col_nb_dofs, false);
358  auxMat.resize(col_nb_dofs, 2, false);
359  int nb_gauss_pts = row_data.getN().size1();
360  for (int gg = 0; gg != nb_gauss_pts; gg++) {
361  // cerr << "gg " << gg << endl;
362  double val = getGaussPts()(2, gg) * getArea();
363  MatrixAdaptor row_diff_n = row_data.getDiffN(gg);
364  MatrixAdaptor col_diff_n = col_data.getDiffN(gg);
365  double an = commonData.aN[gg];
366  ublas::matrix_row<MatrixDouble> grad_u(commonData.gradU, gg);
367  ublas::matrix_row<MatrixDouble> an_pow3_by_grad_u(
369  // cerr << "grad_u " << grad_u << endl;
370  // cerr << "an_pow3_by_grad_u " << an_pow3_by_grad_u << endl;
371  noalias(auxVec) = prod(col_diff_n, an_pow3_by_grad_u);
372  // cerr << "auxVec " << auxVec << endl;
373  noalias(auxMat) = outer_prod(auxVec, grad_u);
374  // cerr << "auxMat " << auxMat << endl;
375  // cerr << row_diff_n << endl;
376  // cerr << col_diff_n << endl;
377  // cerr << prod(row_diff_n,trans(val*(an*col_diff_n+auxMat))) << endl;
378  noalias(kTangent) +=
379  val * prod(row_diff_n, trans(an * col_diff_n - auxMat));
380  // cerr << "kTangent " << kTangent << endl;
381  }
382  CHKERR MatSetValues(getFEMethod()->snes_B, row_nb_dofs,
383  &*row_data.getIndices().data().begin(), col_nb_dofs,
384  &*col_data.getIndices().data().begin(),
385  &*kTangent.data().begin(), ADD_VALUES);
387  }
388  };
389 };
390 } // namespace MinimalSurfaceEquation
391 
392 #endif // __MINIMAL_SURFACE_ELEMENT_HPP__
393 
394 /***************************************************************************/ /**
395  * \defgroup minimal_surface_area Minimal surface area
396  * \ingroup user_modules
397  ******************************************************************************/
EdgeElement feBcEdge
Used to calculate dofs on boundary.
OpGetDataAtGaussPts(const string field_name, CommonData &common_data)
OpAssembleResidaul(const string field_name, CommonData &common_data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Implementation of minimal area element.
SurfaceElement feLhs
To calculate left hand side.
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
VectorDouble normGradU2
size of nb_gauss_pts, norm of gradient
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
MatrixDouble gradU
nb_gauss_pts x 2 gradients at integration pts
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
Definition: Types.hpp:121
OpCalculateCoefficientsAtGaussPts(const string field_name, CommonData &common_data, bool rhs_and_not_lhs)
DataForcesAndSourcesCore::EntData EntData
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:104
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
#define CHKERR
Inline error check.
Definition: definitions.h:596
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Evaluate function values and gradients at Gauss Pts.
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
SurfaceElement feRhs
To calculate right hand side.
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
virtual double evalFunction(double x, double y)
Function on boundary Inherit this class and overload this function to change bc.
Minimal surface equation.
OpAssembleTangent(const string field_name, CommonData &common_data)