v0.14.0
Loading...
Searching...
No Matches
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
78 VectorDouble
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 };
83
84 /** \brief Integrate vector on rhs,
85
86 \f[
87 \mathbf{F} = \int_L \mathbf{N}^\textrm{T} g(x,y) \textrm{d}L
88 \f]
89
90 */
92 : public EdgeElementForcesAndSourcesCore::UserDataOperator {
93 Vec vF;
94 OpAssmebleBcRhs(const string field_name, Vec v_f)
97 vF(v_f) {}
98
99 /** \brief Function on boundary
100 Inherit this class and overload this function to change bc.
101 */
102 virtual double evalFunction(double x, double y) {
103 return sin(2 * M_PI * (x + y));
104 // return sin(M_PI*x);
105 }
106
107 VectorDouble nF;
108 PetscErrorCode doWork(int side, EntityType type,
109 DataForcesAndSourcesCore::EntData &data) {
111
112 int nb_dofs = data.getFieldData().size();
113 if (nb_dofs == 0) {
115 }
116 nF.resize(nb_dofs, false);
117 nF.clear();
118 int nb_gauss_pts = data.getN().size1();
119 for (int gg = 0; gg != nb_gauss_pts; gg++) {
120 double val = getLength() * getGaussPts()(1, gg);
121 double x = getCoordsAtGaussPts()(gg, 0);
122 double y = getCoordsAtGaussPts()(gg, 1);
123 noalias(nF) += val * evalFunction(x, y) * data.getN(gg);
124 }
125 CHKERR VecSetValues(vF, data.getIndices().size(),
126 &*data.getIndices().data().begin(), &nF[0],
127 ADD_VALUES);
129 }
130 };
131
132 /** \brief Integrate vector on lhs,
133
134 \f[
135 \mathbf{A} = \int_L \mathbf{N}^\textrm{T}\mathbf{N} \textrm{d}L
136 \f]
137
138 */
140 : public EdgeElementForcesAndSourcesCore::UserDataOperator {
141 Mat mA;
142 OpAssmebleBcLhs(const string field_name, Mat m_a)
144 field_name, UserDataOperator::OPROWCOL),
145 mA(m_a) {}
146 MatrixDouble K, transK;
147 PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
148 EntityType col_type,
149 DataForcesAndSourcesCore::EntData &row_data,
150 DataForcesAndSourcesCore::EntData &col_data) {
152
153 int row_nb_dofs = row_data.getIndices().size();
154 if (row_nb_dofs == 0) {
156 }
157 int col_nb_dofs = col_data.getIndices().size();
158 if (col_nb_dofs == 0) {
160 }
161 K.resize(row_nb_dofs, col_nb_dofs, false);
162 K.clear();
163 int nb_gauss_pts = row_data.getN().size1();
164 for (int gg = 0; gg != nb_gauss_pts; gg++) {
165 double val = getLength() * getGaussPts()(1, gg);
166 VectorAdaptor row_n = row_data.getN(gg);
167 VectorAdaptor col_n = col_data.getN(gg);
168 noalias(K) += val * outer_prod(row_n, col_n);
169 }
170 CHKERR MatSetValues(mA, row_nb_dofs,
171 &*row_data.getIndices().data().begin(), col_nb_dofs,
172 &*col_data.getIndices().data().begin(),
173 &*K.data().begin(), ADD_VALUES);
174 // Matrix is symmetric, assemble off diagonal part
175 if (row_type != col_type || row_side != col_side) {
176 transK.resize(col_nb_dofs, row_nb_dofs, false);
177 noalias(transK) = trans(K);
178 CHKERR MatSetValues(mA, col_nb_dofs,
179 &*col_data.getIndices().data().begin(), row_nb_dofs,
180 &*row_data.getIndices().data().begin(),
181 &*transK.data().begin(), ADD_VALUES);
182 }
183
185 }
186 };
187
188 /** \brief Evaluate function values and gradients at Gauss Pts
189 */
191 : public FaceElementForcesAndSourcesCore::UserDataOperator {
193 OpGetDataAtGaussPts(const string field_name, CommonData &common_data)
196 commonData(common_data) {}
197
198 PetscErrorCode doWork(int side, EntityType type,
199 DataForcesAndSourcesCore::EntData &data) {
201
202 int nb_dofs = data.getFieldData().size();
203 if (nb_dofs == 0) {
205 }
206 int nb_gauss_pts = data.getN().size1();
207
208 // Element loops over entities, it start from vertices
209 if (type == MBVERTEX) {
210 // Setting size of common data vectors
211 commonData.gradU.resize(nb_gauss_pts, 2);
212 commonData.gradU.clear();
213 }
214
215 for (int gg = 0; gg != nb_gauss_pts; gg++) {
216 MatrixAdaptor grad_shape_fun = data.getDiffN(gg);
217 ublas::matrix_row<MatrixDouble> grad_at_gauss_pt(commonData.gradU, gg);
218 noalias(grad_at_gauss_pt) +=
219 prod(trans(grad_shape_fun), data.getFieldData());
220 }
221
223 }
224 };
225
226 /** \brief Evaluate function values and gradients at Gauss Pts
227 */
229 : public FaceElementForcesAndSourcesCore::UserDataOperator {
233 CommonData &common_data,
234 bool rhs_and_not_lhs)
237 commonData(common_data), rhsAndNotLhs(rhs_and_not_lhs) {}
238
239 PetscErrorCode doWork(int side, EntityType type,
240 DataForcesAndSourcesCore::EntData &data) {
242
243 // Element loops over entities, it start from vertices
244 if (type == MBVERTEX) {
245 int nb_gauss_pts = data.getN().size1();
246 commonData.normGradU2.resize(nb_gauss_pts, false);
247 for (int gg = 0; gg != nb_gauss_pts; gg++) {
248 ublas::matrix_row<MatrixDouble> grad_at_gauss_pt(commonData.gradU,
249 gg);
251 inner_prod(grad_at_gauss_pt, grad_at_gauss_pt);
252 }
253 commonData.aN.resize(nb_gauss_pts, false);
254 for (int gg = 0; gg != nb_gauss_pts; gg++) {
255 commonData.aN[gg] = 1. / sqrt(1 + commonData.normGradU2[gg]);
256 }
257 if (rhsAndNotLhs) {
258 // needed at right hand side when residual is calculated
259 commonData.aNbyGradU.resize(nb_gauss_pts, 2, false);
260 for (int gg = 0; gg != nb_gauss_pts; gg++) {
261 ublas::matrix_row<MatrixDouble> grad_at_gauss_pt(commonData.gradU,
262 gg);
263 for (int rr = 0; rr != 2; rr++) {
264 commonData.aNbyGradU(gg, rr) =
265 commonData.aN[gg] * grad_at_gauss_pt[rr];
266 }
267 }
268 } else {
269 // needed at left hand side when matrix is calculated
270 commonData.aNpow3byGradU.resize(nb_gauss_pts, 2, false);
271 for (int gg = 0; gg != nb_gauss_pts; gg++) {
272 ublas::matrix_row<MatrixDouble> grad_at_gauss_pt(commonData.gradU,
273 gg);
274 for (int rr = 0; rr != 2; rr++) {
275 commonData.aNpow3byGradU(gg, rr) =
276 pow(commonData.aN[gg], 3) * grad_at_gauss_pt[rr];
277 }
278 }
279 }
280 }
281
282 // doVerticesRow = true;
283 // doEdgesRow = false;
284 // doQuadsRow = false;
285 // doTrisRow = false;
286 // doTetsRow = false;
287 // doPrismsRow = false;
288
290 }
291 };
292
293 /** \brief Assemble residual
294 */
296 : public FaceElementForcesAndSourcesCore::UserDataOperator {
298 OpAssembleResidaul(const string field_name, CommonData &common_data)
301 commonData(common_data) {}
302 VectorDouble rEsidual;
303 PetscErrorCode doWork(int side, EntityType type,
304 DataForcesAndSourcesCore::EntData &data) {
306
307 int nb_dofs = data.getIndices().size();
308 if (nb_dofs == 0) {
310 }
311 rEsidual.resize(nb_dofs, false);
312 rEsidual.clear();
313 int nb_gauss_pts = data.getN().size1();
314 for (int gg = 0; gg != nb_gauss_pts; gg++) {
315 double val = getGaussPts()(2, gg) * getArea();
316 ublas::matrix_row<MatrixDouble> an_by_grad_u(commonData.aNbyGradU, gg);
317 noalias(rEsidual) += val * prod(data.getDiffN(gg), an_by_grad_u);
318 }
319 CHKERR VecSetValues(getFEMethod()->snes_f, data.getIndices().size(),
320 &*data.getIndices().data().begin(), &rEsidual[0],
321 ADD_VALUES);
323 }
324 };
325
326 /** \brief Assemble tangent matrix
327 */
329 : public FaceElementForcesAndSourcesCore::UserDataOperator {
331 OpAssembleTangent(const string field_name, CommonData &common_data)
333 field_name, UserDataOperator::OPROWCOL),
334 commonData(common_data) {
335 sYmm = false; // matrix is not symmetric
336 }
337 MatrixDouble kTangent;
338 VectorDouble auxVec;
339 MatrixDouble auxMat;
340 PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
341 EntityType col_type,
342 DataForcesAndSourcesCore::EntData &row_data,
343 DataForcesAndSourcesCore::EntData &col_data) {
345
346 int row_nb_dofs = row_data.getIndices().size();
347 if (row_nb_dofs == 0) {
349 }
350 int col_nb_dofs = col_data.getIndices().size();
351 if (col_nb_dofs == 0) {
353 }
354 kTangent.resize(row_nb_dofs, col_nb_dofs, false);
355 kTangent.clear();
356 auxVec.resize(col_nb_dofs, false);
357 auxMat.resize(col_nb_dofs, 2, false);
358 int nb_gauss_pts = row_data.getN().size1();
359 for (int gg = 0; gg != nb_gauss_pts; gg++) {
360 // cerr << "gg " << gg << endl;
361 double val = getGaussPts()(2, gg) * getArea();
362 MatrixAdaptor row_diff_n = row_data.getDiffN(gg);
363 MatrixAdaptor col_diff_n = col_data.getDiffN(gg);
364 double an = commonData.aN[gg];
365 ublas::matrix_row<MatrixDouble> grad_u(commonData.gradU, gg);
366 ublas::matrix_row<MatrixDouble> an_pow3_by_grad_u(
368 // cerr << "grad_u " << grad_u << endl;
369 // cerr << "an_pow3_by_grad_u " << an_pow3_by_grad_u << endl;
370 noalias(auxVec) = prod(col_diff_n, an_pow3_by_grad_u);
371 // cerr << "auxVec " << auxVec << endl;
372 noalias(auxMat) = outer_prod(auxVec, grad_u);
373 // cerr << "auxMat " << auxMat << endl;
374 // cerr << row_diff_n << endl;
375 // cerr << col_diff_n << endl;
376 // cerr << prod(row_diff_n,trans(val*(an*col_diff_n+auxMat))) << endl;
377 noalias(kTangent) +=
378 val * prod(row_diff_n, trans(an * col_diff_n - auxMat));
379 // cerr << "kTangent " << kTangent << endl;
380 }
381 CHKERR MatSetValues(getFEMethod()->snes_B, row_nb_dofs,
382 &*row_data.getIndices().data().begin(), col_nb_dofs,
383 &*col_data.getIndices().data().begin(),
384 &*kTangent.data().begin(), ADD_VALUES);
386 }
387 };
388};
389} // namespace MinimalSurfaceEquation
390
391#endif // __MINIMAL_SURFACE_ELEMENT_HPP__
392
393/***************************************************************************/ /**
394 * \defgroup minimal_surface_area Minimal surface area
395 * \ingroup user_modules
396 ******************************************************************************/
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr int order
Minimal surface equation.
constexpr auto field_name
MatrixDouble gradU
nb_gauss_pts x 2 gradients at integration pts
VectorDouble normGradU2
size of nb_gauss_pts, norm of gradient
OpAssembleResidaul(const string field_name, CommonData &common_data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpAssembleTangent(const string field_name, CommonData &common_data)
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 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)
virtual double evalFunction(double x, double y)
Function on boundary Inherit this class and overload this function to change bc.
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpCalculateCoefficientsAtGaussPts(const string field_name, CommonData &common_data, bool rhs_and_not_lhs)
Evaluate function values and gradients at Gauss Pts.
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpGetDataAtGaussPts(const string field_name, CommonData &common_data)
Implementation of minimal area element.
EdgeElement feBcEdge
Used to calculate dofs on boundary.
SurfaceElement feRhs
To calculate right hand side.
SurfaceElement feLhs
To calculate left hand side.
Deprecated interface functions.