Calculate Lorentz fore from magnetic field.
Calculate Lorentz fore from magnetic field.It is not attemt to have accurate or realistic model of moving particles in magnetic field. This example was create to show how field evaluator works with precalculated magnetic field.
\[
v_i = (p_{i+1} - p_{i-1}) / (2 \delta t) \\
(p_{i-1} - 2 p_i + p_{i+1}) / \delta t^2 = \frac{q}{m} v_i \times B_i \\
(p_{i-1} - 2 p_i + p_{i+1}) / \delta t^2 = \frac{q}{m} (p_{i+1} - p_{i-1}) /
(2 \delta t) \times B_i \\
(p_{i-1} - 2 p_i + p_{i+1}) / \delta t = \frac{q}{m} (p_{i+1} - p_{i-1})
\times B_i / 2 \\
p_{i+1} / \delta t - p_{i+1} \times B_i / 2= (2 p_i - p_{i-1}) / \delta t -
p_{i-1} \times B_i / 2 \\ p_{i+1} (\mathbf{1} / \delta t - \frac{q}{m}
\mathbf{1} \times B_i / 2 )= (2 p_i - p_{i-1}) / \delta t - \frac{q}{m}
p_{i-1} \times B_i / 2 \]
#include <BasicFiniteElements.hpp>
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
try {
ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
auto moab_comm_wrap =
boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
if (pcomm == NULL)
pcomm =
new ParallelComm(&moab, moab_comm_wrap->get_comm());
PetscBool flg_file;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Lorenz force configure",
"none");
CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"solution.h5m",
ierr = PetscOptionsEnd();
const char *option;
option = "";
DM dm;
CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
CHKERR DMSetType(dm,
"DMMOFEM");
struct MyOpDebug :
public VolOp {
boost::shared_ptr<MatrixDouble>
B;
MyOpDebug(
decltype(
B)
b) :
VolOp(
"MAGNETIC_POTENTIAL", OPROW),
B(
b) {}
if (
type == MBEDGE && side == 0) {
std::cout << "found " << (*B) << endl;
}
}
};
const double dist = 1e-12;
const int nb_random_points = 5000;
const int nb_steps = 10000;
const int mod_step = 10;
const double velocity_scale = 0.1;
const double magnetic_field_scale = 1;
const double scale_box = 0.5;
auto vol_ele = data_at_pts->feMethodPtr.lock();
if (!vol_ele)
"Pointer to element does not exists");
auto get_rule = [&](
int order_row,
int order_col,
int order_data) {
return -1;
};
BoundBox box;
CHKERR data_at_pts->treePtr->get_bounding_box(box);
const double bMin = box.bMin[0];
const double bMax = box.bMax[0];
auto create_vertices = [nb_random_points, bMin, bMax,
auto &arrays_coord) {
ReadUtilIface *iface;
EntityHandle startv;
CHKERR moab.query_interface(iface);
CHKERR iface->get_node_coords(3, nb_random_points, 0, startv,
arrays_coord);
arrays_coord[0], arrays_coord[1], arrays_coord[2]};
verts =
Range(startv, startv + nb_random_points - 1);
for (
int n = 0;
n != nb_random_points; ++
n) {
t_coords(0) = 0;
for (auto ii : {1, 2}) {
t_coords(ii) = scale_box * (bMax - bMin) *
(std::rand() / static_cast<double>(RAND_MAX)) -
bMax * scale_box;
}
++t_coords;
}
};
auto set_positions = [nb_random_points,
dt, velocity_scale](
arrays_coord[0], arrays_coord[1], arrays_coord[2]};
&init_pos(0, 0), &init_pos(1, 0), &init_pos(2, 0)};
for (
int n = 0;
n != nb_random_points; ++
n) {
for (auto ii : {0, 1, 2})
t_velocity(ii) = (rand() / static_cast<double>(RAND_MAX) - 0.5);
t_velocity(
i) /= std::sqrt(t_velocity(
i) * t_velocity(
i));
t_init_coords(
i) = t_coords(
i) +
dt * velocity_scale * t_velocity(
i);
++t_coords;
++t_init_coords;
}
return init_pos;
};
vector<double *> arrays_coord;
CHKERR create_vertices(moab_charged_partices, verts, arrays_coord);
auto init_positions = set_positions(moab, arrays_coord);
auto get_t_coords = [&]() {
arrays_coord[0], arrays_coord[1], arrays_coord[2]);
};
auto get_t_init_coords = [&]() {
&init_positions(0, 0), &init_positions(1, 0), &init_positions(2, 0));
};
auto calc_rhs = [&](auto &t_p, auto &t_init_p, auto &t_B) {
t_rhs(
k) = (2 * t_p(
k) - t_init_p(
k)) /
dt -
return t_rhs;
};
auto calc_lhs = [&](auto &t_B) {
for (auto ii : {0, 1, 2})
return t_lhs;
};
auto is_out = [&](auto &t_p) {
return true;
}
else if (t_p(
i) < bMin) {
return true;
}
return false;
};
auto cache_ptr = boost::make_shared<CacheTuple>();
auto calc_position = [&]() {
auto t_p = get_t_coords();
auto t_init_p = get_t_init_coords();
for (
int n = 0;
n != nb_random_points; ++
n) {
if (is_out(t_p)) {
++t_p;
++t_init_p;
continue;
}
std::array<double, 3> point = {t_p(0), t_p(1), t_p(2)};
data_at_pts->setEvalPoints(point.data(), 1);
point.data(), dist, prb_ptr->
getName(),
"MAGNETIC", data_at_pts,
for (int ii : {0, 1, 2})
t_B(ii) = (*B)(ii, 0);
else
t_B(
i) *= magnetic_field_scale * 0.5;
auto t_rhs = calc_rhs(t_p, t_init_p, t_B);
auto t_lhs = calc_lhs(t_B);
double det;
t_p(
i) = t_inv_lhs(
i,
j) * t_rhs(
j);
++t_p;
++t_init_p;
}
};
for (
int t = 0;
t != nb_steps; ++
t) {
std::string print_step =
"Step : " + boost::lexical_cast<std::string>(
t) +
"\r";
std::cout << print_step << std::flush;
calc_position();
if ((
t % mod_step) == 0) {
CHKERR moab_charged_partices.write_file(
(
"step_" + boost::lexical_cast<std::string>(
t / mod_step) +
".vtk")
.c_str(),
"VTK", "");
}
}
std::cout << endl;
}
return 0;
}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define CATCH_ERRORS
Catch errors.
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ 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< 'n', SPACE_DIM > n
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
FTensor::Index< 'i', SPACE_DIM > i
int main(int argc, char *argv[])
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
DeprecatedCoreInterface Interface
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
constexpr double t
plate stiffness
virtual MoFEMErrorCode cache_problem_entities(const std::string prb_name, CacheTupleWeakPtr cache_ptr)=0
Cache variables.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual int get_comm_rank() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
const VectorDouble & getFieldData() const
get dofs values
Field evaluator interface.
MoFEMErrorCode buildTree3D(boost::shared_ptr< SetPtsData > spd_ptr, const std::string finite_element)
Build spatial tree.
boost::shared_ptr< SPD > getData(const double *ptr=nullptr, const int nb_eval_points=0, const double eps=1e-12, VERBOSITY_LEVELS verb=QUIET)
Get the Data object.
MoFEMErrorCode evalFEAtThePoint3D(const double *const point, const double distance, const std::string problem, const std::string finite_element, boost::shared_ptr< SetPtsData > data_ptr, int lower_rank, int upper_rank, boost::shared_ptr< CacheTuple > cache_ptr, MoFEMTypes bh=MF_EXIST, VERBOSITY_LEVELS verb=QUIET)
Evaluate field at artbitray position.
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
RuleHookFun getRuleHook
Hook to get rule.
Calculate curl of vector field.
keeps basic data about problem
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
VolEle::UserDataOperator VolOp