14 : cOre(const_cast<
Core &>(core)) {}
25 std::vector<RandomFieldData> random_fields,
26 boost::shared_ptr<Range> ents_ptr) {
32 auto get_random_number = [](
auto &range) {
33 const auto r =
static_cast<double>(std::rand()) / RAND_MAX;
34 return range[0] + r * (range[1] - range[0]);
42 VecGetLocalSize(
v, &size);
46 for (
auto &rf : random_fields) {
49 <<
"Set random field " << rf.first <<
" " <<
static_cast<int>(field_id);
52 auto lo = prb_ptr->getNumeredRowDofsPtr()->lower_bound(
54 const auto hi = prb_ptr->getNumeredRowDofsPtr()->upper_bound(
56 for (; lo != hi; ++lo) {
57 const auto idx = (*lo)->getPetscLocalDofIdx();
58 if (idx >= 0 && idx < size) {
59 array[idx] = get_random_number(rf.second);
63 for (
auto pit = ents_ptr->const_pair_begin();
64 pit != ents_ptr->const_pair_end(); ++pit) {
66 prb_ptr->getNumeredRowDofsPtr()->get<
Unique_mi_tag>().lower_bound(
69 prb_ptr->getNumeredRowDofsPtr()->get<
Unique_mi_tag>().upper_bound(
71 for (; lo != hi; ++lo) {
72 const auto idx = (*lo)->getPetscLocalDofIdx();
73 if (idx >= 0 && idx < size) {
74 array[idx] = get_random_number(rf.second);
81 CHKERR VecRestoreArray(
v, &array);
83 CHKERR VecGhostUpdateBegin(
v, INSERT_VALUES, SCATTER_FORWARD);
84 CHKERR VecGhostUpdateEnd(
v, INSERT_VALUES, SCATTER_FORWARD);
91 boost::shared_ptr<FEMethod> pipeline,
99 pipeline->ts_t = time;
100 pipeline->ts_dt = delta_t;
106 auto [
v,
a] =
setPipelineX(pipeline, x, delta_x, delta2_x, delta_t);
114 "Run assemble pipeline");
127 auto get_norm = [](
auto f) {
129 VecNorm(f, NORM_2, &fnorm);
141 boost::shared_ptr<FEMethod> pipeline,
152 pipeline->ts_t = time;
154 pipeline->ts_dt = delta_t;
155 pipeline->ts_a = 1. / delta_t;
156 pipeline->ts_aa = 1. / pow(delta_t, 2);
162 auto [
v,
a] =
setPipelineX(pipeline, x, delta_x, delta2_x, delta_t);
170 "Run assemble pipeline");
174 auto get_norm = [](
auto m) {
176 MatNorm(
m, NORM_1, &fnorm);
193 auto axpy = [&](
auto v0,
auto diff_v,
auto eps) {
195 if (v0.use_count()) {
209 axpy(x, diff_x, -
eps),
211 axpy(delta_x, diff_x, -
eps),
213 axpy(delta2_x, diff_x, -
eps),
215 time, delta_t, cache_ptr);
219 axpy(x, diff_x,
eps),
221 axpy(delta_x, diff_x,
eps),
223 axpy(delta2_x, diff_x,
eps),
225 time, delta_t, cache_ptr);
227 auto calculate_derivative = [
eps](
auto fm1,
auto fp1) {
229 VecGetLocalSize(fp1, &size);
230 VecAXPY(fp1, -1, fm1);
232 VecGetArray(fp1, &array);
233 for (
int j = 0;
j != size; ++
j) {
236 VecRestoreArray(fp1, &array);
240 return calculate_derivative(fm1, fp1);
245 boost::shared_ptr<FEMethod> pipeline_rhs,
252 dm, fe_name, pipeline_rhs, x, delta_x, delta2_x, diff_x, time, delta_t,
255 auto m =
assembleMat(dm, fe_name, pipeline_lhs, x, delta_x, delta2_x, time,
280 VecCopy(delta_x, x_t);
281 VecScale(x_t, 1. / delta_t);
289 pipeline->x_t = PETSC_NULLPTR;
296 VecCopy(delta2_x, x_tt);
297 VecScale(x_tt, 1. / pow(delta_t, 2));
302 pipeline->x_tt = x_tt;
305 pipeline->x_tt = PETSC_NULLPTR;
308 return std::make_pair(x_t, x_tt);
#define MOFEM_TAG_AND_LOG_C(channel, severity, tag, format,...)
Tag and log in channel.
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define CHKERR
Inline error check.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
auto createDMVector(DM dm)
Get smart vector from DM.
auto createDMMatrix(DM dm)
Get smart matrix from DM.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
boost::weak_ptr< CacheTuple > CacheTupleWeakPtr
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto getProblemPtr(DM dm)
get problem pointer from DM
FTensor::Index< 'm', 3 > m
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
Deprecated interface functions.
static UId getLoFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
static UId getHiFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
SmartPetscObj< Mat > assembleMat(SmartPetscObj< DM > dm, std::string fe_name, boost::shared_ptr< FEMethod > pipeline, SmartPetscObj< Vec > x, SmartPetscObj< Vec > delta_x, SmartPetscObj< Vec > delta2_x, double time, double delta_t, CacheTupleWeakPtr cache_ptr)
Assemble the left hand side vector.
SmartPetscObj< Vec > setRandomFields(SmartPetscObj< DM > dm, std::vector< RandomFieldData > random_fields, boost::shared_ptr< Range > ents=nullptr)
Generate random fields.
SmartPetscObj< Vec > checkCentralFiniteDifference(SmartPetscObj< DM > dm, std::string fe_name, boost::shared_ptr< FEMethod > pipeline_rhs, boost::shared_ptr< FEMethod > pipeline_lhs, SmartPetscObj< Vec > x, SmartPetscObj< Vec > delta_x, SmartPetscObj< Vec > delta2_x, SmartPetscObj< Vec > diff_x, double time, double delta_t, double eps, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Check consistency between directional derivative with matrix.
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
OperatorsTester(const MoFEM::Core &core)
SmartPetscObj< Vec > assembleVec(SmartPetscObj< DM > dm, std::string fe_name, boost::shared_ptr< FEMethod > pipeline, SmartPetscObj< Vec > x, SmartPetscObj< Vec > delta_x, SmartPetscObj< Vec > delta2_x, double time, double delta_t, CacheTupleWeakPtr cache_ptr)
Assemble the right hand side vector.
std::pair< SmartPetscObj< Vec >, SmartPetscObj< Vec > > setPipelineX(boost::shared_ptr< FEMethod > pipeline, SmartPetscObj< Vec > x, SmartPetscObj< Vec > delta_x, SmartPetscObj< Vec > delta2_x, double delta_t)
Set vectors, x, x_t, and x_tt to finite element instance.
SmartPetscObj< Vec > directionalCentralFiniteDifference(SmartPetscObj< DM > dm, std::string fe_name, boost::shared_ptr< FEMethod > pipeline, SmartPetscObj< Vec > x, SmartPetscObj< Vec > delta_x, SmartPetscObj< Vec > delta2_x, SmartPetscObj< Vec > diff_x, double time, double delta_t, double eps, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Calculate directional directive using finite difference.
static constexpr Switches CtxSetA
static constexpr Switches CtxSetF
static constexpr Switches CtxSetB
static constexpr Switches CtxSetTime
intrusive_ptr for managing petsc objects
base class for all interface classes