v0.9.0
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 1173 of file MixTransportElement.hpp.

Constructor & Destructor Documentation

◆ OpEvaluateBcOnFluxes()

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

Definition at line 1176 of file MixTransportElement.hpp.

1178  flux_name, UserDataOperator::OPROW),
1179  cTx(ctx) {}
ForcesAndSourcesCore::UserDataOperator UserDataOperator

Member Function Documentation

◆ doWork()

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

Definition at line 1185 of file MixTransportElement.hpp.

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

Member Data Documentation

◆ cTx

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

Definition at line 1175 of file MixTransportElement.hpp.

◆ i

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

Definition at line 1183 of file MixTransportElement.hpp.

◆ Nf

VectorDouble MixTransport::MixTransportElement::OpEvaluateBcOnFluxes::Nf

Definition at line 1182 of file MixTransportElement.hpp.

◆ NN

MatrixDouble MixTransport::MixTransportElement::OpEvaluateBcOnFluxes::NN

Definition at line 1181 of file MixTransportElement.hpp.


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