v0.12.1
Public Member Functions | Public Attributes | List of all members
MixTransport::MixTransportElement::OpEvaluateBcOnFluxes Struct Reference

Evaluate boundary conditions on fluxes. More...

#include <users_modules/tutorials/cor-0to1/src/MixTransportElement.hpp>

Inheritance diagram for MixTransport::MixTransportElement::OpEvaluateBcOnFluxes:
[legend]
Collaboration diagram for MixTransport::MixTransportElement::OpEvaluateBcOnFluxes:
[legend]

Public Member Functions

 OpEvaluateBcOnFluxes (MixTransportElement &ctx, const std::string flux_name)
 
MoFEMErrorCode doWork (int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
 

Public Attributes

MixTransportElementcTx
 
MatrixDouble NN
 
VectorDouble Nf
 
FTensor::Index< 'i', 3 > i
 

Detailed Description

Evaluate boundary conditions on fluxes.

Note that Neumann boundary conditions here are essential. So it is opposite what you find in displacement finite element method.

Here we have to solve for degrees of freedom on boundary such base functions approximate flux.

Definition at line 1126 of file MixTransportElement.hpp.

Constructor & Destructor Documentation

◆ OpEvaluateBcOnFluxes()

MixTransport::MixTransportElement::OpEvaluateBcOnFluxes::OpEvaluateBcOnFluxes ( MixTransportElement ctx,
const std::string  flux_name 
)

Definition at line 1129 of file MixTransportElement.hpp.

1131  flux_name, UserDataOperator::OPROW),
1132  cTx(ctx) {}
ForcesAndSourcesCore::UserDataOperator UserDataOperator

Member Function Documentation

◆ doWork()

MoFEMErrorCode MixTransport::MixTransportElement::OpEvaluateBcOnFluxes::doWork ( int  side,
EntityType  type,
DataForcesAndSourcesCore::EntData data 
)

Definition at line 1138 of file MixTransportElement.hpp.

1139  {
1141  if (data.getFieldData().size() == 0)
1143  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1144  int nb_dofs = data.getFieldData().size();
1145  int nb_gauss_pts = data.getN().size1();
1146  if (3 * nb_dofs != static_cast<int>(data.getN().size2())) {
1147  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1148  "wrong number of dofs");
1149  }
1150  NN.resize(nb_dofs, nb_dofs);
1151  Nf.resize(nb_dofs);
1152  NN.clear();
1153  Nf.clear();
1154 
1155  // Get normal vector. Note that when higher order geometry is set, then
1156  // face element could be curved, i.e. normal can be different at each
1157  // integration point.
1158  double *normal_ptr;
1159  if (getNormalsAtGaussPts().size1() == (unsigned int)nb_gauss_pts) {
1160  // HO geometry
1161  normal_ptr = &getNormalsAtGaussPts(0)[0];
1162  } else {
1163  // Linear geometry, i.e. constant normal on face
1164  normal_ptr = &getNormal()[0];
1165  }
1166  // set tensor from pointer
1167  FTensor::Tensor1<const double *, 3> t_normal(normal_ptr, &normal_ptr[1],
1168  &normal_ptr[2], 3);
1169  // get base functions
1170  auto t_n_hdiv_row = data.getFTensor1N<3>();
1171 
1172  double nrm2 = 0;
1173 
1174  // loop over integration points
1175  for (int gg = 0; gg < nb_gauss_pts; gg++) {
1176 
1177  // get integration point coordinates
1178  const double x = getCoordsAtGaussPts()(gg, 0);
1179  const double y = getCoordsAtGaussPts()(gg, 1);
1180  const double z = getCoordsAtGaussPts()(gg, 2);
1181 
1182  // get flux on fece for given element handle and coordinates
1183  double flux;
1184  CHKERR cTx.getBcOnFluxes(fe_ent, x, y, z, flux);
1185  ;
1186  // get weight for integration rule
1187  double w = getGaussPts()(2, gg);
1188  if (gg == 0) {
1189  nrm2 = sqrt(t_normal(i) * t_normal(i));
1190  }
1191 
1192  // set tensor of rank 0 to matrix NN elements
1193  // loop over base functions on rows and columns
1194  for (int ll = 0; ll != nb_dofs; ll++) {
1195  // get column on shape functions
1197  &data.getVectorN<3>(gg)(0, HVEC0),
1198  &data.getVectorN<3>(gg)(0, HVEC1),
1199  &data.getVectorN<3>(gg)(0, HVEC2), 3);
1200  for (int kk = 0; kk <= ll; kk++) {
1201  NN(ll, kk) += w * t_n_hdiv_row(i) * t_n_hdiv_col(i);
1202  ++t_n_hdiv_col;
1203  }
1204  // right hand side
1205  Nf[ll] += w * t_n_hdiv_row(i) * t_normal(i) * flux / nrm2;
1206  ++t_n_hdiv_row;
1207  }
1208 
1209  // If HO geometry increment t_normal to next integration point
1210  if (getNormalsAtGaussPts().size1() == (unsigned int)nb_gauss_pts) {
1211  ++t_normal;
1212  nrm2 = sqrt(t_normal(i) * t_normal(i));
1213  }
1214  }
1215  // get global dofs indices on element
1216  cTx.bcIndices.insert(data.getIndices().begin(), data.getIndices().end());
1217  // factor matrix
1219  // solve local problem
1220  cholesky_solve(NN, Nf, ublas::lower());
1221 
1222  // set solution to vector
1223  CHKERR VecSetOption(cTx.D0, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
1224  CHKERR VecSetValues(cTx.D0, nb_dofs, &*data.getIndices().begin(),
1225  &*Nf.begin(), INSERT_VALUES);
1226  for (int dd = 0; dd != nb_dofs; ++dd)
1227  data.getFieldDofs()[dd]->getFieldData() = Nf[dd];
1228 
1230  }
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
Definition: cholesky.hpp:221
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
Definition: cholesky.hpp:52
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:440
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:339
@ HVEC0
Definition: definitions.h:179
@ HVEC1
Definition: definitions.h:179
@ HVEC2
Definition: definitions.h:179
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:42
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:409
#define CHKERR
Inline error check.
Definition: definitions.h:528
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
double w(const double g, const double t)
double flux
impulse magnitude
virtual MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x, const double y, const double z, double &flux)
essential (Neumann) boundary condition (set fluxes)

Member Data Documentation

◆ cTx

MixTransportElement& MixTransport::MixTransportElement::OpEvaluateBcOnFluxes::cTx

Definition at line 1128 of file MixTransportElement.hpp.

◆ i

FTensor::Index<'i', 3> MixTransport::MixTransportElement::OpEvaluateBcOnFluxes::i

Definition at line 1136 of file MixTransportElement.hpp.

◆ Nf

VectorDouble MixTransport::MixTransportElement::OpEvaluateBcOnFluxes::Nf

Definition at line 1135 of file MixTransportElement.hpp.

◆ NN

MatrixDouble MixTransport::MixTransportElement::OpEvaluateBcOnFluxes::NN

Definition at line 1134 of file MixTransportElement.hpp.


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