v0.14.0
Loading...
Searching...
No Matches
FieldApproximation.hpp
Go to the documentation of this file.
1/* \file FieldApproximation.hpp
2
3\brief Element to calculate approximation on volume elements
4
5*/
6
7
8
9#ifndef __FILEDAPPROXIMATION_HPP
10#define __FILEDAPPROXIMATION_HPP
11
12using namespace boost::numeric;
13
14/** \brief Finite element for approximating analytical filed on the mesh
15 * \ingroup user_modules
16 */
18
20 const std::string problemName;
21 VolumeElementForcesAndSourcesCore feVolume;
23
25 : mField(m_field), feVolume(m_field), feFace(m_field) {}
26
27 /** \brief Gauss point operators to calculate matrices and vectors
28 *
29 * Function work on volumes (Terahedrons & Bricks)
30 */
31 template <typename FUNEVAL>
34
35 Mat A;
36 std::vector<Vec> &vecF;
38
39 OpApproxVolume(const std::string &field_name, Mat _A,
40 std::vector<Vec> &vec_F, FUNEVAL &function_evaluator)
41 : VolumeElementForcesAndSourcesCore::UserDataOperator(
43 A(_A), vecF(vec_F), functionEvaluator(function_evaluator) {}
44 virtual ~OpApproxVolume() {}
45
46 MatrixDouble NN;
47 MatrixDouble transNN;
48 std::vector<VectorDouble> Nf;
49
50 /** \brief calculate matrix
51 */
52 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
53 EntityType col_type,
54 EntitiesFieldData::EntData &row_data,
55 EntitiesFieldData::EntData &col_data) {
57
58 if (A == PETSC_NULL)
60 if (row_data.getIndices().size() == 0)
62 if (col_data.getIndices().size() == 0)
64
65 const auto &dof_ptr = row_data.getFieldDofs()[0];
66 int rank = dof_ptr->getNbOfCoeffs();
67
68 int nb_row_dofs = row_data.getIndices().size() / rank;
69 int nb_col_dofs = col_data.getIndices().size() / rank;
70
71 NN.resize(nb_row_dofs, nb_col_dofs, false);
72 NN.clear();
73
74 unsigned int nb_gauss_pts = row_data.getN().size1();
75 for (unsigned int gg = 0; gg != nb_gauss_pts; gg++) {
76 double w = getVolume() * getGaussPts()(3, gg);
77 // noalias(NN) += w*outer_prod(row_data.getN(gg),col_data.getN(gg));
78 cblas_dger(CblasRowMajor, nb_row_dofs, nb_col_dofs, w,
79 &row_data.getN()(gg, 0), 1, &col_data.getN()(gg, 0), 1,
80 &*NN.data().begin(), nb_col_dofs);
81 }
82
83 if ((row_type != col_type) || (row_side != col_side)) {
84 transNN.resize(nb_col_dofs, nb_row_dofs, false);
85 ublas::noalias(transNN) = trans(NN);
86 }
87
88 double *data = &*NN.data().begin();
89 double *trans_data = &*transNN.data().begin();
90 VectorInt row_indices, col_indices;
91 row_indices.resize(nb_row_dofs);
92 col_indices.resize(nb_col_dofs);
93
94 for (int rr = 0; rr < rank; rr++) {
95
96 if ((row_data.getIndices().size() % rank) != 0) {
97 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
98 "data inconsistency");
99 }
100
101 if ((col_data.getIndices().size() % rank) != 0) {
102 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
103 "data inconsistency");
104 }
105
106 unsigned int nb_rows;
107 unsigned int nb_cols;
108 int *rows;
109 int *cols;
110
111 if (rank > 1) {
112
113 ublas::noalias(row_indices) = ublas::vector_slice<VectorInt>(
114 row_data.getIndices(),
115 ublas::slice(rr, rank, row_data.getIndices().size() / rank));
116 ublas::noalias(col_indices) = ublas::vector_slice<VectorInt>(
117 col_data.getIndices(),
118 ublas::slice(rr, rank, col_data.getIndices().size() / rank));
119
120 nb_rows = row_indices.size();
121 nb_cols = col_indices.size();
122 rows = &*row_indices.data().begin();
123 cols = &*col_indices.data().begin();
124
125 } else {
126
127 nb_rows = row_data.getIndices().size();
128 nb_cols = col_data.getIndices().size();
129 rows = &*row_data.getIndices().data().begin();
130 cols = &*col_data.getIndices().data().begin();
131 }
132
133 if (nb_rows != NN.size1()) {
134 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
135 "data inconsistency");
136 }
137 if (nb_cols != NN.size2()) {
138 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
139 "data inconsistency");
140 }
141
142 CHKERR MatSetValues(A, nb_rows, rows, nb_cols, cols, data, ADD_VALUES);
143 if ((row_type != col_type) || (row_side != col_side)) {
144 if (nb_rows != transNN.size2()) {
145 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
146 "data inconsistency");
147 }
148 if (nb_cols != transNN.size1()) {
149 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
150 "data inconsistency");
151 }
152 CHKERR MatSetValues(A, nb_cols, cols, nb_rows, rows, trans_data,
153 ADD_VALUES);
154 }
155 }
156
158 }
159
160 /** \brief calculate vector
161 */
162 MoFEMErrorCode doWork(int side, EntityType type,
163 EntitiesFieldData::EntData &data) {
165
166 if (data.getIndices().size() == 0)
168
169 // PetscAttachDebugger();
170
171 const auto &dof_ptr = data.getFieldDofs()[0];
172 unsigned int rank = dof_ptr->getNbOfCoeffs();
173
174 int nb_row_dofs = data.getIndices().size() / rank;
175
176 if (getCoordsAtGaussPts().size1() != data.getN().size1()) {
177 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
178 "data inconsistency");
179 }
180 if (getCoordsAtGaussPts().size2() != 3) {
181 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
182 "data inconsistency");
183 }
184
185 // integration
186 unsigned int nb_gauss_pts = data.getN().size1();
187 for (unsigned int gg = 0; gg != nb_gauss_pts; gg++) {
188
189 const double w = getVolume() * getGaussPts()(3, gg);
190 // intergartion point global positions for linear tetrahedral element
191 const double x = getCoordsAtGaussPts()(gg, 0);
192 const double y = getCoordsAtGaussPts()(gg, 1);
193 const double z = getCoordsAtGaussPts()(gg, 2);
194
195 std::vector<VectorDouble> fun_val;
196
197 fun_val = functionEvaluator(x, y, z);
198
199 if (fun_val.size() != vecF.size()) {
200 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
201 "data inconsistency");
202 }
203
204 Nf.resize(fun_val.size());
205 for (unsigned int lhs = 0; lhs != fun_val.size(); lhs++) {
206
207 if (!gg) {
208 Nf[lhs].resize(data.getIndices().size());
209 Nf[lhs].clear();
210 }
211
212 if (fun_val[lhs].size() != rank) {
213 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
214 "data inconsistency");
215 }
216
217 for (unsigned int rr = 0; rr != rank; rr++) {
218 cblas_daxpy(nb_row_dofs, w * (fun_val[lhs])[rr],
219 &data.getN()(gg, 0), 1, &(Nf[lhs])[rr], rank);
220 }
221 }
222 }
223
224 for (unsigned int lhs = 0; lhs != vecF.size(); lhs++) {
225
226 CHKERR VecSetValues(vecF[lhs], data.getIndices().size(),
227 &data.getIndices()[0], &(Nf[lhs])[0], ADD_VALUES);
228 }
229
231 }
232 };
233
234 /** \brief Gauss point operators to calculate matrices and vectors
235 *
236 * Function work on faces (Triangles & Quads)
237 */
238 template <typename FUNEVAL>
240 : public FaceElementForcesAndSourcesCore::UserDataOperator {
241
242 Mat A;
243 std::vector<Vec> &vecF;
245
246 OpApproxFace(const std::string &field_name, Mat _A, std::vector<Vec> &vec_F,
247 FUNEVAL &function_evaluator)
249 field_name, UserDataOperator::OPROW | UserDataOperator::OPROWCOL),
250 A(_A), vecF(vec_F), functionEvaluator(function_evaluator) {}
251 virtual ~OpApproxFace() {}
252
253 MatrixDouble NN;
254 MatrixDouble transNN;
255 std::vector<VectorDouble> Nf;
256
257 /** \brief calculate matrix
258 */
259 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
260 EntityType col_type,
261 EntitiesFieldData::EntData &row_data,
262 EntitiesFieldData::EntData &col_data) {
264
265 if (A == PETSC_NULL)
267 if (row_data.getIndices().size() == 0)
269 if (col_data.getIndices().size() == 0)
271
272 const auto &dof_ptr = row_data.getFieldDofs()[0];
273 int rank = dof_ptr->getNbOfCoeffs();
274 int nb_row_dofs = row_data.getIndices().size() / rank;
275 int nb_col_dofs = col_data.getIndices().size() / rank;
276 NN.resize(nb_row_dofs, nb_col_dofs, false);
277 NN.clear();
278 unsigned int nb_gauss_pts = row_data.getN().size1();
279 for (unsigned int gg = 0; gg != nb_gauss_pts; gg++) {
280 double w = getGaussPts()(2, gg);
281 if (getNormalsAtGaussPts().size1()) {
282 w *= 0.5 * cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
283 } else {
284 w *= getArea();
285 }
286 // noalias(NN) += w*outer_prod(row_data.getN(gg),col_data.getN(gg));
287 cblas_dger(CblasRowMajor, nb_row_dofs, nb_col_dofs, w,
288 &row_data.getN()(gg, 0), 1, &col_data.getN()(gg, 0), 1,
289 &*NN.data().begin(), nb_col_dofs);
290 }
291
292 if ((row_type != col_type) || (row_side != col_side)) {
293 transNN.resize(nb_col_dofs, nb_row_dofs, false);
294 ublas::noalias(transNN) = trans(NN);
295 }
296
297 double *data = &*NN.data().begin();
298 double *trans_data = &*transNN.data().begin();
299 VectorInt row_indices, col_indices;
300 row_indices.resize(nb_row_dofs);
301 col_indices.resize(nb_col_dofs);
302 for (int rr = 0; rr < rank; rr++) {
303 if ((row_data.getIndices().size() % rank) != 0) {
304 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
305 "data inconsistency");
306 }
307 if ((col_data.getIndices().size() % rank) != 0) {
308 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
309 "data inconsistency");
310 }
311 unsigned int nb_rows;
312 unsigned int nb_cols;
313 int *rows;
314 int *cols;
315 if (rank > 1) {
316 ublas::noalias(row_indices) = ublas::vector_slice<VectorInt>(
317 row_data.getIndices(),
318 ublas::slice(rr, rank, row_data.getIndices().size() / rank));
319 ublas::noalias(col_indices) = ublas::vector_slice<VectorInt>(
320 col_data.getIndices(),
321 ublas::slice(rr, rank, col_data.getIndices().size() / rank));
322 nb_rows = row_indices.size();
323 nb_cols = col_indices.size();
324 rows = &*row_indices.data().begin();
325 cols = &*col_indices.data().begin();
326 } else {
327 nb_rows = row_data.getIndices().size();
328 nb_cols = col_data.getIndices().size();
329 rows = &*row_data.getIndices().data().begin();
330 cols = &*col_data.getIndices().data().begin();
331 }
332 if (nb_rows != NN.size1()) {
333 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
334 "data inconsistency");
335 }
336 if (nb_cols != NN.size2()) {
337 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
338 "data inconsistency");
339 }
340 CHKERR MatSetValues(A, nb_rows, rows, nb_cols, cols, data, ADD_VALUES);
341 if ((row_type != col_type) || (row_side != col_side)) {
342 if (nb_rows != transNN.size2()) {
343 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
344 "data inconsistency");
345 }
346 if (nb_cols != transNN.size1()) {
347 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
348 "data inconsistency");
349 }
350 CHKERR MatSetValues(A, nb_cols, cols, nb_rows, rows, trans_data,
351 ADD_VALUES);
352 }
353 }
355 }
356
357 /** \brief calculate vector
358 */
359 MoFEMErrorCode doWork(int side, EntityType type,
360 EntitiesFieldData::EntData &data) {
362
363 if (data.getIndices().size() == 0)
365
366 const auto &dof_ptr = data.getFieldDofs()[0];
367 unsigned int rank = dof_ptr->getNbOfCoeffs();
368
369 int nb_row_dofs = data.getIndices().size() / rank;
370
371 if (getCoordsAtGaussPts().size1() != data.getN().size1()) {
372 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
373 "data inconsistency");
374 }
375 if (getCoordsAtGaussPts().size2() != 3) {
376 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
377 "data inconsistency");
378 }
379
380 VectorDouble normal(3);
381 VectorDouble tangent1(3);
382 VectorDouble tangent2(3);
383 tangent1.clear();
384 tangent2.clear();
385
386 // integration
387 unsigned int nb_gauss_pts = data.getN().size1();
388 for (unsigned int gg = 0; gg != nb_gauss_pts; gg++) {
389 double w = getGaussPts()(2, gg);
390 w *= 0.5 * cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
391
392 // intergartion point global positions for linear tetrahedral element
393 const double x = getCoordsAtGaussPts()(gg, 0);
394 const double y = getCoordsAtGaussPts()(gg, 1);
395 const double z = getCoordsAtGaussPts()(gg, 2);
396
397 if (getNormalsAtGaussPts().size1()) {
398 noalias(normal) = getNormalsAtGaussPts(gg);
399 for (int dd = 0; dd < 3; dd++) {
400 tangent1[dd] = getTangent1AtGaussPts()(gg, dd);
401 tangent2[dd] = getTangent2AtGaussPts()(gg, dd);
402 }
403 } else {
404 noalias(normal) = getNormal();
405 }
406
407 std::vector<VectorDouble> fun_val;
408 EntityHandle ent = getFEMethod()->numeredEntFiniteElementPtr->getEnt();
409 fun_val = functionEvaluator(ent, x, y, z, normal, tangent1, tangent2);
410 if (fun_val.size() != vecF.size()) {
411 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
412 "data inconsistency");
413 }
414
415 Nf.resize(fun_val.size());
416 for (unsigned int lhs = 0; lhs != fun_val.size(); lhs++) {
417 if (!gg) {
418 Nf[lhs].resize(data.getIndices().size());
419 Nf[lhs].clear();
420 }
421 if (fun_val[lhs].size() != rank) {
422 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
423 "data inconsistency");
424 }
425 for (unsigned int rr = 0; rr != rank; rr++) {
426 cblas_daxpy(nb_row_dofs, w * (fun_val[lhs])[rr],
427 &data.getN()(gg, 0), 1, &(Nf[lhs])[rr], rank);
428 }
429 }
430 }
431
432 for (unsigned int lhs = 0; lhs != vecF.size(); lhs++) {
433
434 CHKERR VecSetValues(vecF[lhs], data.getIndices().size(),
435 &data.getIndices()[0], &(Nf[lhs])[0], ADD_VALUES);
436 }
437
439 }
440 };
441
442 /** \brief Set operators
443 */
444 template <typename FUNEVAL>
445 MoFEMErrorCode setOperatorsVolume(const std::string &field_name, Mat A,
446 std::vector<Vec> &vec_F,
447 FUNEVAL &function_evaluator) {
449 // add operator to calculate F vector
450 feVolume.getOpPtrVector().push_back(
451 new OpApproxVolume<FUNEVAL>(field_name, A, vec_F, function_evaluator));
452 // add operator to calculate A matrix
453 // if(A) {
454 // feVolume.getOpPtrVector().push_back(new
455 // OpApproxVolume<FUNEVAL>(field_name,A,vec_F,function_evaluator));
456 // }
458 }
459
460 /** \brief Set operators
461 */
462 template <typename FUNEVAL>
463 MoFEMErrorCode setOperatorsFace(const std::string &field_name, Mat A,
464 std::vector<Vec> &vec_F,
465 FUNEVAL &function_evaluator) {
467 // add operator to calculate F vector
468 feFace.getOpPtrVector().push_back(
469 new OpApproxFace<FUNEVAL>(field_name, A, vec_F, function_evaluator));
470 // add operator to calculate A matrix
471 // if(A) {
472 // feFace.getOpPtrVector().push_back(new
473 // OpApproxFace<FUNEVAL>(field_name,A,vec_F,function_evaluator));
474 // }
476 }
477
478 /** \brief assemble matrix and vector
479 */
480 template <typename FUNEVAL>
481 MoFEMErrorCode loopMatrixAndVectorVolume(const std::string &problem_name,
482 const std::string &fe_name,
483 const std::string &field_name, Mat A,
484 std::vector<Vec> &vec_F,
485 FUNEVAL &function_evaluator) {
487
488 CHKERR setOperatorsVolume(field_name, A, vec_F, function_evaluator);
489 if (A) {
490 CHKERR MatZeroEntries(A);
491 }
492 // calculate and assemble
493 CHKERR mField.loop_finite_elements(problem_name, fe_name, feVolume);
494 if (A) {
495 CHKERR MatAssemblyBegin(A, MAT_FLUSH_ASSEMBLY);
496 CHKERR MatAssemblyEnd(A, MAT_FLUSH_ASSEMBLY);
497 }
498 for (unsigned int lhs = 0; lhs < vec_F.size(); lhs++) {
499 CHKERR VecAssemblyBegin(vec_F[lhs]);
500 CHKERR VecAssemblyEnd(vec_F[lhs]);
501 }
503 }
504};
505
506#endif //__FILEDAPPROXIMATION_HPP
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
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#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
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
constexpr AssemblyType A
constexpr auto field_name
Gauss point operators to calculate matrices and vectors.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate vector
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
calculate matrix
OpApproxFace(const std::string &field_name, Mat _A, std::vector< Vec > &vec_F, FUNEVAL &function_evaluator)
Gauss point operators to calculate matrices and vectors.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
calculate matrix
OpApproxVolume(const std::string &field_name, Mat _A, std::vector< Vec > &vec_F, FUNEVAL &function_evaluator)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate vector
Finite element for approximating analytical filed on the mesh.
MoFEMErrorCode loopMatrixAndVectorVolume(const std::string &problem_name, const std::string &fe_name, const std::string &field_name, Mat A, std::vector< Vec > &vec_F, FUNEVAL &function_evaluator)
assemble matrix and vector
MoFEM::Interface & mField
VolumeElementForcesAndSourcesCore feVolume
MoFEM::FaceElementForcesAndSourcesCore feFace
MoFEMErrorCode setOperatorsFace(const std::string &field_name, Mat A, std::vector< Vec > &vec_F, FUNEVAL &function_evaluator)
Set operators.
MoFEMErrorCode setOperatorsVolume(const std::string &field_name, Mat A, std::vector< Vec > &vec_F, FUNEVAL &function_evaluator)
Set operators.
FieldApproximationH1(MoFEM::Interface &m_field)
const std::string problemName
Deprecated interface functions.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.