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

Evaluate boundary conditions on fluxes. More...

#include <users_modules/basic_finite_elements/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) {}
ForcesAndSourcesCore::UserDataOperator UserDataOperator

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  }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513
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
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
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
Definition: cholesky.hpp:221
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
#define CHKERR
Inline error check.
Definition: definitions.h:601
virtual MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x, const double y, const double z, double &flux)
essential (Neumann) boundary condition (set fluxes)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
double w(const double g, const double t)
Definition: ContactOps.hpp:160

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: