v0.9.0
Public Member Functions | Public Attributes | List of all members
FieldApproximationH1::OpApproxFace< FUNEVAL > Struct Template Reference

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

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

Inherits UserDataOperator.

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, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
 calculate matrix More...
 
MoFEMErrorCode doWork (int side, EntityType type, DataForcesAndSourcesCore::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 265 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 
)

Definition at line 272 of file FieldApproximation.hpp.

275  field_name, UserDataOperator::OPROW | UserDataOperator::OPROWCOL),
276  A(_A), vecF(vec_F), functionEvaluator(function_evaluator) {}
ForcesAndSourcesCore::UserDataOperator UserDataOperator

◆ ~OpApproxFace()

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

Definition at line 277 of file FieldApproximation.hpp.

277 {}

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,
DataForcesAndSourcesCore::EntData row_data,
DataForcesAndSourcesCore::EntData col_data 
)

calculate matrix

Definition at line 285 of file FieldApproximation.hpp.

288  {
290 
291  if (A == PETSC_NULL)
293  if (row_data.getIndices().size() == 0)
295  if (col_data.getIndices().size() == 0)
297 
298  const auto &dof_ptr = row_data.getFieldDofs()[0];
299  int rank = dof_ptr->getNbOfCoeffs();
300  int nb_row_dofs = row_data.getIndices().size() / rank;
301  int nb_col_dofs = col_data.getIndices().size() / rank;
302  NN.resize(nb_row_dofs, nb_col_dofs, false);
303  NN.clear();
304  unsigned int nb_gauss_pts = row_data.getN().size1();
305  for (unsigned int gg = 0; gg != nb_gauss_pts; gg++) {
306  double w = getGaussPts()(2, gg);
307  if (getNormalsAtGaussPts().size1()) {
308  w *= 0.5 * cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
309  } else {
310  w *= getArea();
311  }
312  // noalias(NN) += w*outer_prod(row_data.getN(gg),col_data.getN(gg));
313  cblas_dger(CblasRowMajor, nb_row_dofs, nb_col_dofs, w,
314  &row_data.getN()(gg, 0), 1, &col_data.getN()(gg, 0), 1,
315  &*NN.data().begin(), nb_col_dofs);
316  }
317 
318  if ((row_type != col_type) || (row_side != col_side)) {
319  transNN.resize(nb_col_dofs, nb_row_dofs, false);
320  ublas::noalias(transNN) = trans(NN);
321  }
322 
323  double *data = &*NN.data().begin();
324  double *trans_data = &*transNN.data().begin();
325  VectorInt row_indices, col_indices;
326  row_indices.resize(nb_row_dofs);
327  col_indices.resize(nb_col_dofs);
328  for (int rr = 0; rr < rank; rr++) {
329  if ((row_data.getIndices().size() % rank) != 0) {
330  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
331  "data inconsistency");
332  }
333  if ((col_data.getIndices().size() % rank) != 0) {
334  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
335  "data inconsistency");
336  }
337  unsigned int nb_rows;
338  unsigned int nb_cols;
339  int *rows;
340  int *cols;
341  if (rank > 1) {
342  ublas::noalias(row_indices) = ublas::vector_slice<VectorInt>(
343  row_data.getIndices(),
344  ublas::slice(rr, rank, row_data.getIndices().size() / rank));
345  ublas::noalias(col_indices) = ublas::vector_slice<VectorInt>(
346  col_data.getIndices(),
347  ublas::slice(rr, rank, col_data.getIndices().size() / rank));
348  nb_rows = row_indices.size();
349  nb_cols = col_indices.size();
350  rows = &*row_indices.data().begin();
351  cols = &*col_indices.data().begin();
352  } else {
353  nb_rows = row_data.getIndices().size();
354  nb_cols = col_data.getIndices().size();
355  rows = &*row_data.getIndices().data().begin();
356  cols = &*col_data.getIndices().data().begin();
357  }
358  if (nb_rows != NN.size1()) {
359  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
360  "data inconsistency");
361  }
362  if (nb_cols != NN.size2()) {
363  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
364  "data inconsistency");
365  }
366  CHKERR MatSetValues(A, nb_rows, rows, nb_cols, cols, data, ADD_VALUES);
367  if ((row_type != col_type) || (row_side != col_side)) {
368  if (nb_rows != transNN.size2()) {
369  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
370  "data inconsistency");
371  }
372  if (nb_cols != transNN.size1()) {
373  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
374  "data inconsistency");
375  }
376  CHKERR MatSetValues(A, nb_cols, cols, nb_rows, rows, trans_data,
377  ADD_VALUES);
378  }
379  }
381  }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
void cblas_dger(const enum CBLAS_ORDER order, const int M, const int N, const double alpha, const double *X, const int incX, const double *Y, const int incY, double *A, const int lda)
Definition: cblas_dger.c:12
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
double cblas_dnrm2(const int N, const double *X, const int incX)
Definition: cblas_dnrm2.c:12
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
#define CHKERR
Inline error check.
Definition: definitions.h:596
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:71
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ doWork() [2/2]

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

calculate vector

Definition at line 385 of file FieldApproximation.hpp.

386  {
388 
389  if (data.getIndices().size() == 0)
391 
392  const auto &dof_ptr = data.getFieldDofs()[0];
393  unsigned int rank = dof_ptr->getNbOfCoeffs();
394 
395  int nb_row_dofs = data.getIndices().size() / rank;
396 
397  if (getCoordsAtGaussPts().size1() != data.getN().size1()) {
398  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
399  "data inconsistency");
400  }
401  if (getCoordsAtGaussPts().size2() != 3) {
402  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
403  "data inconsistency");
404  }
405 
406  VectorDouble normal(3);
407  VectorDouble tangent1(3);
408  VectorDouble tangent2(3);
409  tangent1.clear();
410  tangent2.clear();
411 
412  // integration
413  unsigned int nb_gauss_pts = data.getN().size1();
414  for (unsigned int gg = 0; gg != nb_gauss_pts; gg++) {
415  double x, y, z, w;
416  w = getGaussPts()(2, gg);
417  if (getNormalsAtGaussPts().size1()) {
418  w *= 0.5 * cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
419  } else {
420  w *= getArea();
421  }
422 
423  if (getHoCoordsAtGaussPts().size1() == nb_gauss_pts) {
424  // intergation points global positions if higher order geometry is
425  // given
426  x = getHoCoordsAtGaussPts()(gg, 0);
427  y = getHoCoordsAtGaussPts()(gg, 1);
428  z = getHoCoordsAtGaussPts()(gg, 2);
429  } else {
430  // intergartion point global positions for linear tetrahedral element
431  x = getCoordsAtGaussPts()(gg, 0);
432  y = getCoordsAtGaussPts()(gg, 1);
433  z = getCoordsAtGaussPts()(gg, 2);
434  }
435 
436  if (getNormalsAtGaussPts().size1()) {
437  noalias(normal) = getNormalsAtGaussPts(gg);
438  for (int dd = 0; dd < 3; dd++) {
439  tangent1[dd] = getTangent1AtGaussPts()(gg, dd);
440  tangent2[dd] = getTangent2AtGaussPts()(gg, dd);
441  }
442  } else {
443  noalias(normal) = getNormal();
444  }
445 
446  std::vector<VectorDouble> fun_val;
447  EntityHandle ent = getFEMethod()->numeredEntFiniteElementPtr->getEnt();
448  fun_val = functionEvaluator(ent, x, y, z, normal, tangent1, tangent2);
449  if (fun_val.size() != vecF.size()) {
450  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
451  "data inconsistency");
452  }
453 
454  Nf.resize(fun_val.size());
455  for (unsigned int lhs = 0; lhs != fun_val.size(); lhs++) {
456  if (!gg) {
457  Nf[lhs].resize(data.getIndices().size());
458  Nf[lhs].clear();
459  }
460  if (fun_val[lhs].size() != rank) {
461  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
462  "data inconsistency");
463  }
464  for (unsigned int rr = 0; rr != rank; rr++) {
465  cblas_daxpy(nb_row_dofs, w * (fun_val[lhs])[rr],
466  &data.getN()(gg, 0), 1, &(Nf[lhs])[rr], rank);
467  }
468  }
469  }
470 
471  for (unsigned int lhs = 0; lhs != vecF.size(); lhs++) {
472 
473  CHKERR VecSetValues(vecF[lhs], data.getIndices().size(),
474  &data.getIndices()[0], &(Nf[lhs])[0], ADD_VALUES);
475  }
476 
478  }
void cblas_daxpy(const int N, const double alpha, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_daxpy.c:11
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
double cblas_dnrm2(const int N, const double *X, const int incX)
Definition: cblas_dnrm2.c:12
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
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72

Member Data Documentation

◆ A

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

Definition at line 268 of file FieldApproximation.hpp.

◆ functionEvaluator

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

Definition at line 270 of file FieldApproximation.hpp.

◆ Nf

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

Definition at line 281 of file FieldApproximation.hpp.

◆ NN

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

Definition at line 279 of file FieldApproximation.hpp.

◆ transNN

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

Definition at line 280 of file FieldApproximation.hpp.

◆ vecF

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

Definition at line 269 of file FieldApproximation.hpp.


The documentation for this struct was generated from the following file: