v0.13.2
Loading...
Searching...
No Matches
Public Member Functions | Public Attributes | List of all members
FieldApproximationH1::OpApproxFace< FUNEVAL > Struct Template Reference
User modules

Gauss point operators to calculate matrices and vectors. More...

#include <users_modules/basic_finite_elements/src/FieldApproximation.hpp>

Inheritance diagram for FieldApproximationH1::OpApproxFace< FUNEVAL >:
[legend]
Collaboration diagram for FieldApproximationH1::OpApproxFace< FUNEVAL >:
[legend]

Public Member Functions

 OpApproxFace (const std::string &field_name, Mat _A, std::vector< Vec > &vec_F, FUNEVAL &function_evaluator)
 
virtual ~OpApproxFace ()
 
MoFEMErrorCode doWork (int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
 calculate matrix More...
 
MoFEMErrorCode doWork (int side, EntityType type, EntitiesFieldData::EntData &data)
 calculate vector More...
 

Public Attributes

Mat A
 
std::vector< Vec > & vecF
 
FUNEVAL & functionEvaluator
 
MatrixDouble NN
 
MatrixDouble transNN
 
std::vector< VectorDouble > Nf
 

Detailed Description

template<typename FUNEVAL>
struct FieldApproximationH1::OpApproxFace< FUNEVAL >

Gauss point operators to calculate matrices and vectors.

Function work on faces (Triangles & Quads)

Definition at line 239 of file FieldApproximation.hpp.

Constructor & Destructor Documentation

◆ OpApproxFace()

template<typename FUNEVAL >
FieldApproximationH1::OpApproxFace< FUNEVAL >::OpApproxFace ( const std::string &  field_name,
Mat  _A,
std::vector< Vec > &  vec_F,
FUNEVAL &  function_evaluator 
)
inline

Definition at line 246 of file FieldApproximation.hpp.

248 : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
249 field_name, UserDataOperator::OPROW | UserDataOperator::OPROWCOL),
250 A(_A), vecF(vec_F), functionEvaluator(function_evaluator) {}
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13

◆ ~OpApproxFace()

template<typename FUNEVAL >
virtual FieldApproximationH1::OpApproxFace< FUNEVAL >::~OpApproxFace ( )
inlinevirtual

Definition at line 251 of file FieldApproximation.hpp.

251{}

Member Function Documentation

◆ doWork() [1/2]

template<typename FUNEVAL >
MoFEMErrorCode FieldApproximationH1::OpApproxFace< FUNEVAL >::doWork ( int  row_side,
int  col_side,
EntityType  row_type,
EntityType  col_type,
EntitiesFieldData::EntData &  row_data,
EntitiesFieldData::EntData &  col_data 
)
inline

calculate matrix

Definition at line 259 of file FieldApproximation.hpp.

262 {
263 MoFEMFunctionBegin;
264
265 if (A == PETSC_NULL)
266 MoFEMFunctionReturnHot(0);
267 if (row_data.getIndices().size() == 0)
268 MoFEMFunctionReturnHot(0);
269 if (col_data.getIndices().size() == 0)
270 MoFEMFunctionReturnHot(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 }
354 MoFEMFunctionReturn(0);
355 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEMFunctionBegin
#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
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::Types::VectorInt
UBlasVector< int > VectorInt
Definition: Types.hpp:67
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1631
OpContactTools::w
double w(const double g, const double t)
Definition: ContactOperators.hpp:109

◆ doWork() [2/2]

template<typename FUNEVAL >
MoFEMErrorCode FieldApproximationH1::OpApproxFace< FUNEVAL >::doWork ( int  side,
EntityType  type,
EntitiesFieldData::EntData &  data 
)
inline

calculate vector

Definition at line 359 of file FieldApproximation.hpp.

360 {
361 MoFEMFunctionBegin;
362
363 if (data.getIndices().size() == 0)
364 MoFEMFunctionReturnHot(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
438 MoFEMFunctionReturn(0);
439 }
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1576

Member Data Documentation

◆ A

template<typename FUNEVAL >
Mat FieldApproximationH1::OpApproxFace< FUNEVAL >::A

Definition at line 242 of file FieldApproximation.hpp.

◆ functionEvaluator

template<typename FUNEVAL >
FUNEVAL& FieldApproximationH1::OpApproxFace< FUNEVAL >::functionEvaluator

Definition at line 244 of file FieldApproximation.hpp.

◆ Nf

template<typename FUNEVAL >
std::vector<VectorDouble> FieldApproximationH1::OpApproxFace< FUNEVAL >::Nf

Definition at line 255 of file FieldApproximation.hpp.

◆ NN

template<typename FUNEVAL >
MatrixDouble FieldApproximationH1::OpApproxFace< FUNEVAL >::NN

Definition at line 253 of file FieldApproximation.hpp.

◆ transNN

template<typename FUNEVAL >
MatrixDouble FieldApproximationH1::OpApproxFace< FUNEVAL >::transNN

Definition at line 254 of file FieldApproximation.hpp.

◆ vecF

template<typename FUNEVAL >
std::vector<Vec>& FieldApproximationH1::OpApproxFace< FUNEVAL >::vecF

Definition at line 243 of file FieldApproximation.hpp.

The documentation for this struct was generated from the following file:
Generated by Doxygen 1.9.5 and hosted at University of Glasgow