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

Evaluate boundary conditions on fluxes. More...

#include <users_modules/tutorials/other_tutorials/mix_transport/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 1171 of file MixTransportElement.hpp.

Constructor & Destructor Documentation

◆ OpEvaluateBcOnFluxes()

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

Definition at line 1174 of file MixTransportElement.hpp.

1176  flux_name, UserDataOperator::OPROW),
1177  cTx(ctx) {}

Member Function Documentation

◆ doWork()

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

Definition at line 1183 of file MixTransportElement.hpp.

1184  {
1186  if (data.getFieldData().size() == 0)
1188  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1189  int nb_dofs = data.getFieldData().size();
1190  int nb_gauss_pts = data.getN().size1();
1191  if (3 * nb_dofs != static_cast<int>(data.getN().size2())) {
1192  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1193  "wrong number of dofs");
1194  }
1195  NN.resize(nb_dofs, nb_dofs);
1196  Nf.resize(nb_dofs);
1197  NN.clear();
1198  Nf.clear();
1199 
1200  // Get normal vector. Note that when higher order geometry is set, then
1201  // face element could be curved, i.e. normal can be different at each
1202  // integration point.
1203  double *normal_ptr;
1204  if (getNormalsAtGaussPts().size1() == (unsigned int)nb_gauss_pts) {
1205  // HO geometry
1206  normal_ptr = &getNormalsAtGaussPts(0)[0];
1207  } else {
1208  // Linear geometry, i.e. constant normal on face
1209  normal_ptr = &getNormal()[0];
1210  }
1211  // set tensor from pointer
1212  FTensor::Tensor1<const double *, 3> t_normal(normal_ptr, &normal_ptr[1],
1213  &normal_ptr[2], 3);
1214  // get base functions
1215  auto t_n_hdiv_row = data.getFTensor1N<3>();
1216 
1217  double nrm2 = 0;
1218 
1219  // loop over integration points
1220  for (int gg = 0; gg < nb_gauss_pts; gg++) {
1221 
1222  // get integration point coordinates
1223  double x, y, z;
1224  if (getNormalsAtGaussPts().size1() == (unsigned int)nb_gauss_pts) {
1225  x = getHoCoordsAtGaussPts()(gg, 0);
1226  y = getHoCoordsAtGaussPts()(gg, 1);
1227  z = getHoCoordsAtGaussPts()(gg, 2);
1228  } else {
1229  x = getCoordsAtGaussPts()(gg, 0);
1230  y = getCoordsAtGaussPts()(gg, 1);
1231  z = getCoordsAtGaussPts()(gg, 2);
1232  }
1233 
1234  // get flux on fece for given element handle and coordinates
1235  double flux;
1236  CHKERR cTx.getBcOnFluxes(fe_ent, x, y, z, flux);
1237  ;
1238  // get weight for integration rule
1239  double w = getGaussPts()(2, gg);
1240  if (gg == 0) {
1241  nrm2 = sqrt(t_normal(i) * t_normal(i));
1242  }
1243 
1244  // set tensor of rank 0 to matrix NN elements
1245  // loop over base functions on rows and columns
1246  for (int ll = 0; ll != nb_dofs; ll++) {
1247  // get column on shape functions
1249  &data.getVectorN<3>(gg)(0, HVEC0),
1250  &data.getVectorN<3>(gg)(0, HVEC1),
1251  &data.getVectorN<3>(gg)(0, HVEC2), 3);
1252  for (int kk = 0; kk <= ll; kk++) {
1253  NN(ll, kk) += w * t_n_hdiv_row(i) * t_n_hdiv_col(i);
1254  ++t_n_hdiv_col;
1255  }
1256  // right hand side
1257  Nf[ll] += w * t_n_hdiv_row(i) * t_normal(i) * flux / nrm2;
1258  ++t_n_hdiv_row;
1259  }
1260 
1261  // If HO geometry increment t_normal to next integration point
1262  if (getNormalsAtGaussPts().size1() == (unsigned int)nb_gauss_pts) {
1263  ++t_normal;
1264  nrm2 = sqrt(t_normal(i) * t_normal(i));
1265  }
1266  }
1267  // get global dofs indices on element
1268  cTx.bcIndices.insert(data.getIndices().begin(), data.getIndices().end());
1269  // factor matrix
1271  // solve local problem
1272  cholesky_solve(NN, Nf, ublas::lower());
1273 
1274  // set solution to vector
1275  CHKERR VecSetOption(cTx.D0, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
1276  CHKERR VecSetValues(cTx.D0, nb_dofs, &*data.getIndices().begin(),
1277  &*Nf.begin(), INSERT_VALUES);
1278  for (int dd = 0; dd != nb_dofs; ++dd)
1279  data.getFieldDofs()[dd]->getFieldData() = Nf[dd];
1280 
1282  }

Member Data Documentation

◆ cTx

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

Definition at line 1173 of file MixTransportElement.hpp.

◆ i

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

Definition at line 1181 of file MixTransportElement.hpp.

◆ Nf

VectorDouble MixTransport::MixTransportElement::OpEvaluateBcOnFluxes::Nf

Definition at line 1180 of file MixTransportElement.hpp.

◆ NN

MatrixDouble MixTransport::MixTransportElement::OpEvaluateBcOnFluxes::NN

Definition at line 1179 of file MixTransportElement.hpp.


The documentation for this struct was generated from the following file:
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
EntityHandle
MixTransport::MixTransportElement::OpEvaluateBcOnFluxes::Nf
VectorDouble Nf
Definition: MixTransportElement.hpp:1180
MixTransport::MixTransportElement::OpEvaluateBcOnFluxes::i
FTensor::Index< 'i', 3 > i
Definition: MixTransportElement.hpp:1181
FTensor::Tensor1
Definition: Tensor1_value.hpp:9
OpContactTools::w
double w(const double g, const double t)
Definition: ContactOperators.hpp:141
HVEC1
@ HVEC1
Definition: definitions.h:255
MixTransport::MixTransportElement::bcIndices
set< int > bcIndices
Definition: MixTransportElement.hpp:95
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:604
MixTransport::MixTransportElement::getBcOnFluxes
virtual MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x, const double y, const double z, double &flux)
essential (Neumann) boundary condition (set fluxes)
Definition: MixTransportElement.hpp:178
MixTransport::MixTransportElement::OpEvaluateBcOnFluxes::cTx
MixTransportElement & cTx
Definition: MixTransportElement.hpp:1173
MoFEM::FaceElementForcesAndSourcesCoreSwitch::UserDataOperator
FaceElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
Definition: FaceElementForcesAndSourcesCore.hpp:327
MixTransport::MixTransportElement::D0
Vec D0
Definition: MixTransportElement.hpp:482
MixTransport::MixTransportElement::OpEvaluateBcOnFluxes::NN
MatrixDouble NN
Definition: MixTransportElement.hpp:1179
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:123
cholesky_decompose
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
cholesky_solve
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
Definition: cholesky.hpp:221
Poisson2DTraditionalDDOperators::VecSetValues
constexpr auto VecSetValues
Definition: poisson_dd_H.hpp:19
HVEC2
@ HVEC2
Definition: definitions.h:255
convert.int
int
Definition: convert.py:66
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
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
HVEC0
@ HVEC0
Definition: definitions.h:255