2133 {
2135
2137
2139 int nb_integration_pts = getGaussPts().size2();
2140 int nb_base_functions = data.
getN().size2();
2141
2142 double time = getFEMethod()->ts_t;
2145 }
2146
2147#ifndef NDEBUG
2148 if (this->locF.size() != nb_dofs)
2150 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
2151#endif
2152
2153 auto integrate_rhs = [&](auto &bc, auto calc_tau, double time_scale) {
2155
2156 auto t_val = getFTensor1FromPtr<3>(&*bc.vals.begin());
2158 auto t_w = getFTensor0IntegrationWeight();
2159 auto t_coords = getFTensor1CoordsAtGaussPts();
2160 auto t_normal = getFTensor1NormalsAtGaussPts();
2161
2163
2164 for (int gg = 0; gg != nb_integration_pts; ++gg) {
2165
2166 double a = sqrt(t_normal(
i) * t_normal(
i));
2168 const auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
2169 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
2170 int rr = 0;
2171 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
2173 (time_scale *
a * t_w * t_row_base * tau) * (t_val(
i) *
scale);
2174 ++t_row_base;
2175 ++t_f;
2176 }
2177
2178 for (; rr != nb_base_functions; ++rr)
2179 ++t_row_base;
2180 ++t_w;
2181 ++t_coords;
2182 ++t_normal;
2183 }
2185 };
2186
2187
2189 for (
auto &bc : *(
bcData)) {
2190 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2191
2192 double time_scale = 1;
2195 }
2196
2198 if (nb_dofs) {
2199
2200 if (std::regex_match(bc.blockName, std::regex(".*COOK.*"))) {
2202 y -= 44;
2203 y /= (60 - 44);
2204 return -y * (y - 1) / 0.25;
2205 };
2206 CHKERR integrate_rhs(bc, calc_tau, time_scale);
2207 } else {
2209 bc, [](double, double, double) { return 1; }, time_scale);
2210 }
2211 }
2212 }
2213 }
2215}
#define FTENSOR_INDEX(DIM, I)
#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.
FTensor::Index< 'i', SPACE_DIM > i
static PetscBool physicalTimeFlg
static double currentPhysicalTime
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
Get DOF values on entity.