1#ifndef __ELECPHYSOPERATORS_HPP__
2#define __ELECPHYSOPERATORS_HPP__
16using EntData = DataForcesAndSourcesCore::EntData;
29const double gma = 0.002;
32const double mu1 = 0.20;
33const double mu2 = 0.30;
52 jac.resize(3, 3,
false);
64 int nb_dofs = data.getIndices().size();
72 int size2 = data.getN().size2();
73 if (3 * nb_dofs !=
static_cast<int>(data.getN().size2()))
75 "wrong number of dofs");
76 nN.resize(nb_dofs, nb_dofs,
false);
77 nF.resize(nb_dofs,
false);
81 auto t_row_tau = data.getFTensor1N<3>();
98 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
100 nrm2 = sqrt(t_normal(
i) * t_normal(
i));
102 const double a = t_w * vol;
103 for (
int rr = 0; rr != nb_dofs; rr++) {
105 &data.getVectorN<3>(gg)(0,
HVEC0),
106 &data.getVectorN<3>(gg)(0,
HVEC1),
107 &data.getVectorN<3>(gg)(0,
HVEC2), 3);
109 for (
int cc = 0; cc != nb_dofs; cc++) {
110 nN(rr, cc) +=
a * (t_row_tau(
i) * t_normal(
i)) *
111 (t_col_tau(
i) * t_normal(
i)) / (nrm2 * nrm2);
119 nrm2 = sqrt(t_normal(
i) * t_normal(
i));
127 for (
auto &dof : data.getFieldDofs()) {
128 dof->getFieldData() =
nF[dof->getEntDofIdx()];
158 int nb_dofs = data.getFieldData().size();
164 if (nb_dofs !=
static_cast<int>(data.getN().size2()))
166 "wrong number of dofs");
167 nN.resize(nb_dofs, nb_dofs,
false);
168 nF.resize(nb_dofs,
false);
172 auto t_row_mass = data.getFTensor0N();
176 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
177 const double a = t_w * vol;
179 for (
int rr = 0; rr != nb_dofs; rr++) {
180 auto t_col_mass = data.getFTensor0N(gg, 0);
182 for (
int cc = 0; cc != nb_dofs; cc++) {
183 nN(rr, cc) +=
a * t_row_mass * t_col_mass;
194 for (
auto &dof : data.getFieldDofs()) {
195 dof->getFieldData() =
nF[dof->getEntDofIdx()];
208 typedef boost::function<
double(
const double,
const double,
const double)>
211 boost::shared_ptr<PreviousData> &data_u,
212 boost::shared_ptr<PreviousData> &data_v,
219 boost::shared_ptr<PreviousData>
dataU;
220 boost::shared_ptr<PreviousData>
dataV;
228 int nb_dofs = data.getFieldData().size();
231 if (nb_dofs !=
static_cast<int>(data.getN().size2()))
233 "wrong number of dofs");
234 nN.resize(nb_dofs, nb_dofs,
false);
235 nF.resize(nb_dofs,
false);
245 auto t_row_mass = data.getFTensor0N();
249 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
250 const double a = t_w * vol;
252 for (
int rr = 0; rr != nb_dofs; rr++) {
253 auto t_col_mass = data.getFTensor0N(gg, 0);
254 nF[rr] +=
a * vn * t_row_mass;
255 for (
int cc = 0; cc != nb_dofs; cc++) {
256 nN(rr, cc) +=
a * t_row_mass * t_col_mass;
269 for (
auto &dof : data.getFieldDofs()) {
270 dof->getFieldData() =
nF[dof->getEntDofIdx()];
287 typedef boost::function<
double(
const double,
const double)>
290 boost::shared_ptr<PreviousData> &data_u,
291 boost::shared_ptr<PreviousData> &data_v,
306 const int nb_dofs = data.getIndices().size();
308 if (nb_dofs !=
static_cast<int>(data.getN().size2()))
310 "wrong number of dofs");
311 vecF.resize(nb_dofs,
false);
312 mat.resize(nb_dofs, nb_dofs,
false);
315 const int nb_integration_pts =
getGaussPts().size2();
320 auto t_row_v_base = data.getFTensor0N();
325 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
327 double rhs =
rhsU(t_val_u, t_val_v);
328 const double a = vol * t_w;
330 for (
int rr = 0; rr != nb_dofs; ++rr) {
331 auto t_col_v_base = data.getFTensor0N(gg, 0);
332 vecF[rr] +=
a * rhs * t_row_v_base;
333 for (
int cc = 0; cc != nb_dofs; ++cc) {
334 mat(rr, cc) +=
a * t_row_v_base * t_col_v_base;
356 boost::shared_ptr<PreviousData>
dataU;
357 boost::shared_ptr<PreviousData>
dataV;
420 boost::shared_ptr<PreviousData> &data)
426 const int nb_dofs = data.getIndices().size();
430 vecF.resize(nb_dofs,
false);
433 const int nb_integration_pts =
getGaussPts().size2();
434 auto t_flux_value = getFTensor1FromMat<3>(
commonData->flux_values);
436 auto t_tau_base = data.getFTensor1N<3>();
438 auto t_tau_grad = data.getFTensor2DiffN<3, 3>();
443 for (
int gg = 0; gg < nb_integration_pts; ++gg) {
445 const double K_inv = 1. /
D_tilde;
446 const double a = vol * t_w;
447 for (
int rr = 0; rr < nb_dofs; ++rr) {
449 t_tau_grad(0, 0) + t_tau_grad(1, 1) + t_tau_grad(2, 2);
450 vecF[rr] += (K_inv * t_tau_base(
i) * t_flux_value(
i) -
451 div_base * t_mass_value) *
477 boost::shared_ptr<PreviousData> &data_u,
487 const int nb_dofs = data.getIndices().size();
491 vecF.resize(nb_dofs,
false);
493 const int nb_integration_pts =
getGaussPts().size2();
498 auto t_row_v_base = data.getFTensor0N();
509 double duration = 5.0;
511 if (
T -
dt < c_time && c_time <=
T + duration) {
520 for (
int gg = 0; gg < nb_integration_pts; ++gg) {
521 const double a = vol * t_w;
522 for (
int rr = 0; rr < nb_dofs; ++rr) {
523 vecF[rr] += t_row_v_base * (t_dot_u + t_div_u - stim) *
a;
539 boost::shared_ptr<PreviousData>
dataU;
540 boost::shared_ptr<PreviousData>
dataV;
559 const int nb_row_dofs = row_data.getIndices().size();
560 const int nb_col_dofs = col_data.getIndices().size();
562 if (nb_row_dofs && nb_col_dofs) {
564 mat.resize(nb_row_dofs, nb_col_dofs,
false);
566 const int nb_integration_pts =
getGaussPts().size2();
568 auto t_row_tau_base = row_data.getFTensor1N<3>();
573 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
574 const double a = vol * t_w;
575 const double K_inv = 1. /
D_tilde;
576 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
577 auto t_col_tau_base = col_data.getFTensor1N<3>(gg, 0);
578 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
579 mat(rr, cc) += (K_inv * t_row_tau_base(
i) * t_col_tau_base(
i)) *
a;
588 if (row_side != col_side || row_type != col_type) {
589 transMat.resize(nb_col_dofs, nb_row_dofs,
false);
615 const int nb_row_dofs = row_data.getIndices().size();
616 const int nb_col_dofs = col_data.getIndices().size();
618 if (nb_row_dofs && nb_col_dofs) {
620 mat.resize(nb_row_dofs, nb_col_dofs,
false);
622 const int nb_integration_pts =
getGaussPts().size2();
624 auto t_row_tau_base = row_data.getFTensor1N<3>();
626 auto t_row_tau_grad = row_data.getFTensor2DiffN<3, 3>();
629 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
630 const double a = vol * t_w;
632 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
633 auto t_col_v_base = col_data.getFTensor0N(gg, 0);
634 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
635 double div_row_base = t_row_tau_grad(0, 0) + t_row_tau_grad(1, 1) +
636 t_row_tau_grad(2, 2);
637 mat(rr, cc) += -(div_row_base * t_col_v_base) *
a;
667 const int nb_row_dofs = row_data.getIndices().size();
668 const int nb_col_dofs = col_data.getIndices().size();
670 if (nb_row_dofs && nb_col_dofs) {
671 mat.resize(nb_row_dofs, nb_col_dofs,
false);
673 const int nb_integration_pts =
getGaussPts().size2();
675 auto t_row_v_base = row_data.getFTensor0N();
678 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
679 const double a = vol * t_w;
680 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
681 auto t_col_tau_grad = col_data.getFTensor2DiffN<3, 3>(gg, 0);
682 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
683 double div_col_base = t_col_tau_grad(0, 0) + t_col_tau_grad(1, 1) + t_col_tau_grad(2, 2);
684 mat(rr, cc) += (t_row_v_base * div_col_base) *
a;
714 const int nb_row_dofs = row_data.getIndices().size();
715 const int nb_col_dofs = col_data.getIndices().size();
716 if (nb_row_dofs && nb_col_dofs) {
718 mat.resize(nb_row_dofs, nb_col_dofs,
false);
720 const int nb_integration_pts =
getGaussPts().size2();
722 auto t_row_v_base = row_data.getFTensor0N();
728 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
729 const double a = vol * t_w;
731 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
732 auto t_col_v_base = col_data.getFTensor0N(gg, 0);
733 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
734 mat(rr, cc) += (ts_a * t_row_v_base * t_col_v_base) *
a;
744 if (row_side != col_side || row_type != col_type) {
745 transMat.resize(nb_col_dofs, nb_row_dofs,
false);
854 boost::shared_ptr<PostProcVolumeOnRefinedMesh> &post_proc)
863 "out_level_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
872 boost::shared_ptr<PostProcVolumeOnRefinedMesh>
postProc;
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
const int save_every_nth_step
DataForcesAndSourcesCore::EntData EntData
FTensor::Index< 'i', 3 > i
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Monitor(MPI_Comm &comm, const int &rank, SmartPetscObj< DM > &dm, boost::shared_ptr< PostProcVolumeOnRefinedMesh > &post_proc)
boost::shared_ptr< PostProcVolumeOnRefinedMesh > postProc
MoFEMErrorCode operator()()
function is run for every finite element
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode preProcess()
function is run at the beginning of loop
OpAssembleLhsTauTau(std::string flux_field)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpAssembleLhsTauV(std::string flux_field, std::string mass_field)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpAssembleLhsVTau(std::string mass_field, std::string flux_field)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpAssembleLhsVV(std::string mass_field)
OpAssembleSlowRhsV(std::string mass_field, boost::shared_ptr< PreviousData > &data_u, boost::shared_ptr< PreviousData > &data_v, Feval_u rhs_u)
boost::function< double(const double, const double)> Feval_u
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< PreviousData > dataU
boost::shared_ptr< PreviousData > dataV
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< PreviousData > commonData
OpAssembleStiffRhsTau(std::string flux_field, boost::shared_ptr< PreviousData > &data)
boost::shared_ptr< PreviousData > dataV
boost::shared_ptr< PreviousData > dataU
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpAssembleStiffRhsV(std::string mass_field, boost::shared_ptr< PreviousData > &data_u, Range &stim_region)
Range & essential_bd_ents
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpEssentialBC(std::string flux_field, Range &essential_bd_ents)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpInitialMass(const std::string &mass_field, Range &inner_surface, double init_val)
OpSolveRecovery(const std::string &mass_field, boost::shared_ptr< PreviousData > &data_u, boost::shared_ptr< PreviousData > &data_v, Method runge_kutta4)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::function< double(const double, const double, const double)> Method
boost::shared_ptr< PreviousData > dataV
boost::shared_ptr< PreviousData > dataU
bool sYmm
If true assume that matrix is symmetric structure.
structure for User Loop Methods on finite elements
default operator for TRI element
MatrixDouble & getNormalsAtGaussPts()
if higher order geometry return normals at Gauss pts.
VectorDouble & getNormal()
get triangle normal
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
intrusive_ptr for managing petsc objects
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
PetscInt ts_step
time step number
Volume finite element base.