24#ifndef __MINIMAL_SURFACE_ELEMENT_HPP__
25#define __MINIMAL_SURFACE_ELEMENT_HPP__
92 :
public EdgeElementForcesAndSourcesCore::UserDataOperator {
103 return sin(2 * M_PI * (x + y));
109 DataForcesAndSourcesCore::EntData &data) {
112 int nb_dofs = data.getFieldData().size();
116 nF.resize(nb_dofs,
false);
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);
125 CHKERR VecSetValues(
vF, data.getIndices().size(),
126 &*data.getIndices().data().begin(), &
nF[0],
140 :
public EdgeElementForcesAndSourcesCore::UserDataOperator {
149 DataForcesAndSourcesCore::EntData &row_data,
150 DataForcesAndSourcesCore::EntData &col_data) {
153 int row_nb_dofs = row_data.getIndices().size();
154 if (row_nb_dofs == 0) {
157 int col_nb_dofs = col_data.getIndices().size();
158 if (col_nb_dofs == 0) {
161 K.resize(row_nb_dofs, col_nb_dofs,
false);
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);
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);
175 if (row_type != col_type || row_side != col_side) {
176 transK.resize(col_nb_dofs, row_nb_dofs,
false);
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);
191 :
public FaceElementForcesAndSourcesCore::UserDataOperator {
199 DataForcesAndSourcesCore::EntData &data) {
202 int nb_dofs = data.getFieldData().size();
206 int nb_gauss_pts = data.getN().size1();
209 if (type == MBVERTEX) {
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());
229 :
public FaceElementForcesAndSourcesCore::UserDataOperator {
234 bool rhs_and_not_lhs)
240 DataForcesAndSourcesCore::EntData &data) {
244 if (type == MBVERTEX) {
245 int nb_gauss_pts = data.getN().size1();
247 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
251 inner_prod(grad_at_gauss_pt, grad_at_gauss_pt);
254 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
260 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
263 for (
int rr = 0; rr != 2; rr++) {
271 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
274 for (
int rr = 0; rr != 2; rr++) {
296 :
public FaceElementForcesAndSourcesCore::UserDataOperator {
304 DataForcesAndSourcesCore::EntData &data) {
307 int nb_dofs = data.getIndices().size();
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();
317 noalias(
rEsidual) += val * prod(data.getDiffN(gg), an_by_grad_u);
319 CHKERR VecSetValues(getFEMethod()->snes_f, data.getIndices().size(),
320 &*data.getIndices().data().begin(), &
rEsidual[0],
329 :
public FaceElementForcesAndSourcesCore::UserDataOperator {
342 DataForcesAndSourcesCore::EntData &row_data,
343 DataForcesAndSourcesCore::EntData &col_data) {
346 int row_nb_dofs = row_data.getIndices().size();
347 if (row_nb_dofs == 0) {
350 int col_nb_dofs = col_data.getIndices().size();
351 if (col_nb_dofs == 0) {
354 kTangent.resize(row_nb_dofs, col_nb_dofs,
false);
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++) {
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);
366 ublas::matrix_row<MatrixDouble> an_pow3_by_grad_u(
370 noalias(
auxVec) = prod(col_diff_n, an_pow3_by_grad_u);
378 val * prod(row_diff_n, trans(an * col_diff_n -
auxMat));
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);
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
Minimal surface equation.
constexpr auto field_name
Keep date between operators.
MatrixDouble gradU
nb_gauss_pts x 2 gradients at integration pts
VectorDouble normGradU2
size of nb_gauss_pts, norm of gradient
VectorDouble aN
size of nb_gauss_pts,
MatrixDouble aNpow3byGradU
nb_gauss_pts x 2,
MatrixDouble aNbyGradU
nb_gauss_pts x 2,
int getRule(int order)
Set integrate rule.
EdgeElement(MoFEM::Interface &m_field)
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)
Integrate vector on lhs,.
OpAssmebleBcLhs(const string field_name, Mat m_a)
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Integrate vector on rhs,.
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.
OpAssmebleBcRhs(const string field_name, Vec v_f)
Evaluate function values and gradients at Gauss Pts.
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)
SurfaceElement(MoFEM::Interface &m_field)
int getRule(int order)
Set integrate rule.
Implementation of minimal area element.
MoFEM::Interface & mField
EdgeElement feBcEdge
Used to calculate dofs on boundary.
SurfaceElement feRhs
To calculate right hand side.
SurfaceElement feLhs
To calculate left hand side.
MinimalSurfaceElement(MoFEM::Interface &m_field)
Deprecated interface functions.