v0.9.1
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)
Deprecated interface functions.
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:75
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
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:513
MatrixDouble gradU
nb_gauss_pts x 2 gradients at integration pts
ForcesAndSourcesCore::UserDataOperator UserDataOperator
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
Definition: Types.hpp:125
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:108
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
#define CHKERR
Inline error check.
Definition: definitions.h:601
constexpr int order
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:412
SurfaceElement feRhs
To calculate right hand side.
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:73
virtual double evalFunction(double x, double y)
Function on boundary Inherit this class and overload this function to change bc.
Minimal surface equation.
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
OpAssembleTangent(const string field_name, CommonData &common_data)