59 template <
int spacedim,
class Value>
66 "GMSH mesh with data. Can be different from actual computational mesh.")
68 "Section where to find the field.\n Some sections are specific to file format: "
69 "point_data/node_data, cell_data/element_data, -/element_node_data, native/-.\n"
70 "If not given by a user, we try to find the field in all sections, but we report an error "
71 "if it is found in more than one section.")
73 "The values of the Field are read from the ```$ElementData``` section with field name given by this key.")
75 "Default value is set on elements which values have not been listed in the mesh data file.")
77 "Definition of the unit of all times defined in the mesh data file.")
79 "This key allows reading field data from the mesh data file shifted in time. Considering the time 't', field descriptor with time 'T', "
80 "time shift 'S', then if 't > T', we read the time frame 't + S'.")
82 IT::Default(
"\"equivalent_mesh\""),
"Type of interpolation applied to the input spatial data.\n"
83 "The default value 'equivalent_mesh' assumes the data being constant on elements living on the same mesh "
84 "as the computational mesh, but possibly with different numbering. In the case of the same numbering, "
85 "the user can set 'identical_mesh' to omit algorithm for guessing node and element renumbering. "
86 "Alternatively, in case of different input mesh, several interpolation algorithms are available.")
88 "Distinguishes bulk / boundary FieldFE.")
92 template <
int spacedim,
class Value>
96 "Specify the section in mesh input file where field data is listed.\nSome sections are specific to file format.")
98 .
add_value(OutputTime::DiscreteSpace::ELEM_DATA,
"element_data",
"cell_data (VTK) / element_data (GMSH)")
100 .
add_value(OutputTime::DiscreteSpace::NATIVE_DATA,
"native_data",
"native_data (only for VTK)")
104 template <
int spacedim,
class Value>
107 return it::Selection(
"interpolation",
"Specify interpolation of the input data from its input mesh to the computational mesh.")
108 .
add_value(DataInterpolation::identic_msh,
"identic_mesh",
"Topology and indices of nodes and elements of"
109 "the input mesh and the computational mesh are identical. "
110 "This interpolation is typically used for GMSH input files containing only the field values without "
111 "explicit mesh specification.")
112 .
add_value(DataInterpolation::equivalent_msh,
"equivalent_mesh",
"Topologies of the input mesh and the computational mesh "
113 "are the same, the node and element numbering may differ. "
114 "This interpolation can be used also for VTK input data.")
115 .
add_value(DataInterpolation::gauss_p0,
"P0_gauss",
"Topologies of the input mesh and the computational mesh may differ. "
116 "Constant values on the elements of the computational mesh are evaluated using the Gaussian quadrature of the fixed order 4, "
117 "where the quadrature points and their values are found in the input mesh and input data using the BIH tree search."
119 .
add_value(DataInterpolation::interp_p0,
"P0_intersection",
"Topologies of the input mesh and the computational mesh may differ. "
120 "Can be applied only for boundary fields. For every (boundary) element of the computational mesh the "
121 "intersection with the input mesh is computed. Constant values on the elements of the computational mesh "
122 "are evaluated as the weighted average of the (constant) values on the intersecting elements of the input mesh.")
126 template <
int spacedim,
class Value>
128 Input::register_class< FieldFE<spacedim, Value>,
unsigned int >(
"FieldFE") +
133 template <
int spacedim,
class Value>
137 boundary_domain_(false), fe_values_(4)
143 template<
int spacedim,
class Value>
148 unsigned int position = 0;
150 if (multifield_arr.
size() > 1)
151 while (
index_ != position) {
156 if (field_rec.
val<std::string>(
"TYPE") ==
"FieldFE") {
159 if (discretization == OutputTime::DiscreteSpace::NATIVE_DATA) {
160 std::shared_ptr< FieldFE<spacedim, Value> > field_fe = std::make_shared< FieldFE<spacedim, Value> >(field.
n_comp());
162 field_fe->init_from_input( *
it, init_data );
174 template <
int spacedim,
class Value>
178 if (dof_values.
size()==0) {
186 this->fill_fe_item<0>();
187 this->fill_fe_item<1>();
188 this->fill_fe_item<2>();
189 this->fill_fe_item<3>();
190 this->
fe_ =
dh_->ds()->fe();
192 this->fill_fe_system_data<0>(block_index);
193 this->fill_fe_system_data<1>(block_index);
194 this->fill_fe_system_data<2>(block_index);
195 this->fill_fe_system_data<3>(block_index);
214 template <
int spacedim,
class Value>
228 unsigned int last_element_idx = -1;
231 unsigned int range_bgn=0, range_end=0;
238 for (
unsigned int i_data = reg_chunk_begin; i_data < reg_chunk_end; ++i_data) {
240 if (elm_idx != last_element_idx) {
243 cell =
dh_->cell_accessor_from_element( elm_idx );
245 last_element_idx = elm_idx;
246 range_bgn = this->
fe_item_[elm.
dim()].range_begin_;
253 for (
unsigned int i_dof=range_bgn, i_cdof=0; i_dof<range_end; i_dof++, i_cdof++) {
256 data_cache.
set(i_data) = mat_value;
261 template <
int spacedim,
class Value>
264 std::shared_ptr<EvalPoints> eval_points = cache_map.
eval_points();
265 std::array<Quadrature, 4> quads{
QGauss(0, 1), this->init_quad<1>(eval_points), this->init_quad<2>(eval_points), this->init_quad<3>(eval_points)};
273 template <
int spacedim,
class Value>
274 template <
unsigned int dim>
278 for (
unsigned int k=0; k<eval_points->size(dim); k++)
279 quad.
set(k) = eval_points->local_point<dim>(k);
284 template <
int spacedim,
class Value>
308 template <
int spacedim,
class Value>
312 ASSERT(
field_name_ !=
"").error(
"Uninitialized FieldFE, did you call init_from_input()?\n");
323 WarningOut().fmt(
"Source mesh of FieldFE '{}' is not compatible with target mesh.\nInterpolation of input data will be changed to 'P0_gauss'.\n",
326 }
else if (this->
interpolation_ == DataInterpolation::interp_p0) {
329 WarningOut().fmt(
"Interpolation 'P0_intersection' of FieldFE '{}' can't be used on bulk region.\nIt will be changed to 'P0_gauss'.\n",
334 if (
dh_ ==
nullptr) {
344 template <
int spacedim,
class Value>
349 switch (this->
value_.n_rows() * this->value_.n_cols()) {
368 std::shared_ptr<DOFHandlerMultiDim> dh_par = std::make_shared<DOFHandlerMultiDim>(
const_cast<MeshBase &
>(*mesh) );
369 std::shared_ptr<DiscreteSpace> ds = std::make_shared<EqualOrderDiscreteSpace>( &
const_cast<MeshBase &
>(*mesh), fe);
370 dh_par->distribute_dofs(ds);
373 this->fill_fe_item<0>();
374 this->fill_fe_item<1>();
375 this->fill_fe_item<2>();
376 this->fill_fe_item<3>();
377 this->
fe_ =
dh_->ds()->fe();
384 template <
int spacedim,
class Value>
388 ASSERT(
field_name_ !=
"").error(
"Uninitialized FieldFE, did you call init_from_input()?\n");
392 unsigned int n_components = this->
value_.n_rows() * this->
value_.n_cols();
395 double read_time = (time.
end()+time_shift) / time_unit_coef;
397 unsigned int n_entities, bdr_shift;
398 bool is_native = (this->
discretization_ == OutputTime::DiscreteSpace::NATIVE_DATA);
400 n_components *=
dh_->max_elem_dofs();
407 n_entities = reader_mesh->n_elements() + reader_mesh->bc_mesh()->n_elements();
408 bdr_shift = reader_mesh->n_elements();
413 auto header = reader->find_header(header_query);
415 header, n_entities, n_components, bdr_shift);
419 }
else if (this->
interpolation_==DataInterpolation::identic_msh) {
421 }
else if (this->
interpolation_==DataInterpolation::equivalent_msh) {
435 template <
int spacedim,
class Value>
438 static const unsigned int quadrature_order = 4;
443 unsigned int quadrature_size=0;
445 unsigned int elem_count;
451 QGauss quad(3, quadrature_order);
452 q_points.resize(quad.
size());
453 q_weights.resize(quad.
size());
456 for (
auto cell :
dh_->own_range()) {
457 auto ele = cell.elm();
458 std::fill(elem_value.begin(), elem_value.end(), 0.0);
459 switch (cell.dim()) {
462 q_points[0] = *ele.node(0);
466 quadrature_size = compute_fe_quadrature<1>(q_points, q_weights, ele, quadrature_order);
469 quadrature_size = compute_fe_quadrature<2>(q_points, q_weights, ele, quadrature_order);
472 quadrature_size = compute_fe_quadrature<3>(q_points, q_weights, ele, quadrature_order);
475 searched_elements.clear();
476 source_mesh->get_bih_tree().find_bounding_box(ele.bounding_box(), searched_elements);
478 auto r_idx = cell.elm().region_idx().idx();
479 std::string reg_name = cell.elm().region().label();
480 for (
unsigned int i=0; i<quadrature_size; ++i) {
481 std::fill(sum_val.begin(), sum_val.end(), 0.0);
486 switch (elm->
dim()) {
488 contains = arma::norm(*elm.
node(0) - q_points[i], 2) < 4*std::numeric_limits<double>::epsilon();
504 unsigned int index = sum_val.size() * (*it);
505 for (
unsigned int j=0; j < sum_val.size(); j++) {
512 if (elem_count > 0) {
513 for (
unsigned int j=0; j < sum_val.size(); j++) {
514 elem_value[j] += (sum_val[j] / elem_count) * q_weights[i];
520 loc_dofs = cell.get_loc_dof_indices();
522 ASSERT_LE(loc_dofs.n_elem, elem_value.size());
523 for (
unsigned int i=0; i < elem_value.size(); i++) {
531 template <
int spacedim,
class Value>
537 double total_measure;
540 for (
auto elm :
dh_->mesh()->elements_range()) {
541 if (elm.dim() == 3) {
542 THROW( ExcInvalidElemeDim() << EI_ElemIdx(elm.idx()) );
545 double epsilon = 4* numeric_limits<double>::epsilon() * elm.measure();
546 auto r_idx = elm.region_idx().idx();
547 std::string reg_name = elm.region().label();
550 if (elm.dim() == 0) {
551 searched_elements.clear();
552 source_mesh->get_bih_tree().find_point(*elm.node(0), searched_elements);
555 searched_elements.clear();
556 source_mesh->get_bih_tree().find_bounding_box(bb, searched_elements);
570 if (source_elm->
dim() == 3) {
574 arma::vec::fixed<3> real_point = *elm.node(0);
578 measure = (std::fabs(arma::sum( unit_point )-1) <= 1e-14
579 && arma::min( unit_point ) >= 0)
606 if (measure > epsilon) {
607 unsigned int index =
value.size() * (*it);
608 for (
unsigned int i=0; i <
value.size(); i++) {
611 total_measure += measure;
617 if (total_measure > epsilon) {
622 for (
unsigned int i=0; i <
value.size(); i++) {
626 WarningOut().fmt(
"Processed element with idx {} is out of source mesh!\n", elm.idx());
634 template <
int spacedim,
class Value>
638 unsigned int data_vec_i, vec_inc;
655 for (
auto cell :
dh_->own_range()) {
656 LocDofVec loc_dofs = cell.get_loc_dof_indices();
658 int source_idx = source_target_vec[cell.elm_idx()];
661 data_vec_i = source_idx;
664 data_vec_i = (source_idx + shift) *
dh_->max_elem_dofs();
667 auto r_idx = cell.elm().region_idx().idx();
668 std::string reg_name = cell.elm().region().label();
669 for (
unsigned int i=0; i<loc_dofs.n_elem; ++i) {
672 ++count_vector[ loc_dofs[i] ];
673 data_vec_i += vec_inc;
756 template <
int spacedim,
class Value>
760 double loc_values[n_vals];
763 for (
auto dh_cell :
dh_->own_range()) {
764 LocDofVec loc_dofs = dh_cell.get_loc_dof_indices();
765 for (i=0; i<loc_dofs.n_elem; ++i) loc_values[i] =
data_vec_.
get( loc_dofs[i] );
766 for ( ; i<n_vals; ++i) loc_values[i] = numeric_limits<double>::signaling_NaN();
767 output_data_cache.
store_value( dh_cell.local_idx(), loc_values );
775 template <
int spacedim,
class Value>
782 template <
int spacedim,
class Value>
789 template <
int spacedim,
class Value>
796 template <
int spacedim,
class Value>
805 return_val = (*input_data_cache_)[i_cache_el];
808 actual_compute_region_error =
RegionValueErr(region_name, elm_idx, return_val);
832 template <
int spacedim,
class Value>
#define ASSERT_PERMANENT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
#define ASSERT_LT(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
#define ASSERT_GT(a, b)
Definition of comparative assert macro (Greater Than) only for debug mode.
#define ASSERT_PTR(ptr)
Definition of assert macro checking non-null pointer (PTR) only for debug mode.
#define ASSERT_LE(a, b)
Definition of comparative assert macro (Less or Equal) only for debug mode.
ArrayMatSet set(uint index)
Class represents boundary part of mesh.
Bounding box in 3d ambient space.
unsigned int compute(std::vector< IPAux > &IP13s)
Computes intersection points for 1D-3D intersection. Computes lower dimensional CIs abscissa vs tetra...
void init()
Initializes lower dimensional objects. Sets correctly the pointers to Plucker coordinates and product...
Class for 2D-2D intersections.
void compute(IntersectionAux< 2, 3 > &intersection)
Computes intersection points for 2D-3D intersection. Computes lower dimensional CIs: 1) 3x triangle s...
void init()
Initializes lower dimensional objects. Sets correctly the pointers to Plucker coordinates and product...
Cell accessor allow iterate over DOF handler cells.
LocDofVec get_loc_dof_indices() const
Returns the local indices of dofs associated to the cell on the local process.
unsigned int dim() const
Return dimension of element appropriate to cell.
NodeAccessor< 3 > node(unsigned int ni) const
unsigned int idx() const
We need this method after replacing Region by RegionIdx, and movinf RegionDB instance into particular...
Directing class of FieldValueCache.
unsigned int region_chunk_end(unsigned int region_patch_idx) const
Return end position of region chunk in FieldValueCache.
std::shared_ptr< EvalPoints > eval_points() const
Getter of eval_points object.
unsigned int region_chunk_begin(unsigned int region_patch_idx) const
Return begin position of region chunk in FieldValueCache.
unsigned int region_idx_from_chunk_position(unsigned int chunk_pos) const
Return begin position of region chunk specified by position in map.
const EvalPointData & eval_point_data(unsigned int point_idx) const
Return item of eval_point_data_ specified by its position.
unsigned int n_comp() const
void set_dof_handler_hash(std::size_t hash)
unsigned int n_dofs_per_element() const
void store_value(unsigned int idx, const T *value)
Compound finite element on dim dimensional simplex.
const std::vector< std::shared_ptr< FiniteElement< dim > > > & fe() const
double unit_conversion_coefficient_
Coeficient of conversion of user-defined unit.
bool is_constant_in_space_
Flag detects that field is only dependent on time.
void init_unit_conversion_coefficient(const Input::Record &rec, const struct FieldAlgoBaseInitData &init_data)
Init value of unit_conversion_coefficient_ from input.
Value value_
Last value, prevents passing large values (vectors) by value.
Common abstract parent of all Field<...> classes.
const std::string & input_name() const
std::pair< double, double > limits() const
FieldFlag::Flags get_flags() const
FieldCommon & units(const UnitSI &units)
Set basic units of the field.
unsigned int n_comp() const
Field< spacedim, Value >::FieldBasePtr create_field(Input::Record rec, const FieldCommon &field) override
std::shared_ptr< DOFHandlerMultiDim > conc_dof_handler_
void cache_reinit(const ElementCacheMap &cache_map) override
Armor::ArmaMat< typename Value::element_type, Value::NRows_, Value::NCols_ > handle_fe_shape(unsigned int dim, unsigned int i_dof, unsigned int i_qp)
double get_scaled_value(int i_cache_el, unsigned int elm_idx, const std::string ®ion_name, RegionValueErr &actual_compute_region_error)
MixedPtr< FiniteElement > fe_
VectorMPI data_vec_
Store data of Field.
static const Input::Type::Selection & get_interp_selection_input_type()
unsigned int data_size() const
std::vector< RegionValueErr > region_value_err_
Set holds data of valid / invalid element values on all regions.
void local_to_ghost_data_scatter_begin()
Call begin scatter functions (local to ghost) on data vector.
void calculate_element_values()
Calculate data of equivalent_mesh interpolation or native data on input over all elements of target m...
virtual void init_from_input(const Input::Record &rec, const struct FieldAlgoBaseInitData &init_data)
virtual ~FieldFE()
Destructor.
void local_to_ghost_data_scatter_end()
Call end scatter functions (local to ghost) on data vector.
static const Input::Type::Record & get_input_type()
Implementation.
std::vector< FEValues< spacedim > > fe_values_
List of FEValues objects of dimensions 0,1,2,3 used for value calculation.
void interpolate_gauss()
Interpolate data (use Gaussian distribution) over all elements of target mesh.
void make_dof_handler(const MeshBase *mesh)
Create DofHandler object.
Quadrature init_quad(std::shared_ptr< EvalPoints > eval_points)
Initialize FEValues object of given dimension.
OutputTime::DiscreteSpace discretization_
Specify section where to find the field data in input mesh file.
FilePath reader_file_
mesh reader file
std::array< FEItem, 4 > fe_item_
Holds specific data of field evaluation over all dimensions.
bool set_time(const TimeStep &time) override
void interpolate_intersection()
Interpolate data (use intersection library) over all elements of target mesh.
VectorMPI set_fe_data(std::shared_ptr< DOFHandlerMultiDim > dh, VectorMPI dof_values=VectorMPI::sequential(0), unsigned int block_index=FieldFE< spacedim, Value >::undef_uint)
double default_value_
Default value of element if not set in mesh data file.
void cache_update(FieldValueCache< typename Value::element_type > &data_cache, ElementCacheMap &cache_map, unsigned int region_patch_idx) override
bool boundary_domain_
Is set in set_mesh method. Value true means, that we accept only boundary element accessors in the va...
DataInterpolation interpolation_
Specify type of FE data interpolation.
Input::Record in_rec_
Accessor to Input::Record.
FieldFlag::Flags flags_
Field flags.
void native_data_to_cache(ElementDataCache< double > &output_data_cache)
std::string field_name_
field name read from input
FieldFE(unsigned int n_comp=0)
std::shared_ptr< DOFHandlerMultiDim > dh_
DOF handler object.
std::shared_ptr< EquivalentMeshMap > source_target_mesh_elm_map_
Maps element indices from computational mesh to the source (data).
static const Input::Type::Selection & get_disc_selection_input_type()
ElementDataCache< double >::CacheData input_data_cache_
Input ElementDataCache is stored in set_time and used in all evaluation and interpolation methods.
void set_mesh(const Mesh *mesh) override
static constexpr Mask declare_input
The field can be set from input. The key in input field descriptor is declared. (default on)
static constexpr Mask equation_input
The field is data parameter of the owning equation. (default on)
std::shared_ptr< FieldBaseType > FieldBasePtr
Dedicated class for storing path to input and output files.
constexpr bool match(Mask mask) const
Internal auxiliary class representing intersection object of simplex<dimA> and simplex<dimB>.
Class represents intersection of two elements.
double compute_measure() const override
Computes the relative measure of intersection object.
static bool contains_point(arma::vec point, ElementAccessor< 3 > elm)
Test if element contains given point.
static BaryPoint project_real_to_unit(const RealPoint &point, const ElementMap &map)
static ElementMap element_map(ElementAccessor< 3 > elm)
Base class for Mesh and BCMesh.
const RegionDB & region_db() const
static const unsigned int undef_idx
BCMesh * bc_mesh() const override
Implement MeshBase::bc_mesh(), getter of boundary mesh.
The class for outputting data during time.
Symmetric Gauss-Legendre quadrature formulae on simplices.
Base class for quadrature rules on simplices in arbitrary dimensions.
unsigned int size() const
Returns number of quadrature points.
Armor::Array< double >::ArrayMatSet set(uint i)
static std::shared_ptr< EquivalentMeshMap > identic_mesh_map(const FilePath &file_path, Mesh *computational_mesh)
static std::shared_ptr< EquivalentMeshMap > eqivalent_mesh_map(const FilePath &file_path, Mesh *computational_mesh)
static std::shared_ptr< Mesh > get_mesh(const FilePath &file_path)
static std::shared_ptr< BaseMeshReader > get_reader(const FilePath &file_path)
unsigned int size() const
static const Input::Type::Tuple & get_input_time_type(double lower_bound=-std::numeric_limits< double >::max(), double upper_bound=std::numeric_limits< double >::max())
Representation of one time step..
double read_time(Input::Iterator< Input::Tuple > time_it, double default_time=std::numeric_limits< double >::quiet_NaN()) const
double read_coef(Input::Iterator< Input::Record > unit_it) const
static Input::Type::Default get_input_default()
static const Input::Type::Record & get_input_type()
double get(unsigned int pos) const
Return value on given position.
void normalize(unsigned int pos, double divisor)
Normalize value on given position.
void set(unsigned int pos, double val)
Set value on given position.
void local_to_ghost_end()
local_to_ghost_{begin,end} updates the ghost values on neighbouring processors from local values
void add(unsigned int pos, double val)
Add value to item on given position.
static VectorMPI sequential(unsigned int size)
void local_to_ghost_begin()
local_to_ghost_{begin,end} updates the ghost values on neighbouring processors from local values
unsigned int size() const
Return size of output data.
Fundamental simplicial intersections.
Definitions of basic Lagrangean finite elements with polynomial shape functions.
Class FESystem for compound finite elements.
MixedPtr< FESystem > mixed_fe_system(MixedPtr< FiniteElement > fe, Args &&... args)
#define INSTANCE_ALL(field)
#define FLOW123D_FORCE_LINK_IN_CHILD(x)
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
arma::Col< IntIdx > LocDofVec
int LongIdx
Define type that represents indices of large arrays (elements, nodes, dofs etc.)
Internal class representing intersection object.
Classes with algorithms for computation of intersections of meshes.
static constexpr bool value
#define WarningOut()
Macro defining 'warning' record of log.
Class MappingP1 implements the affine transformation of the unit cell onto the actual cell.
typename arma::Mat< Type >::template fixed< nr, nc > ArmaMat
Definitions of particular quadrature rules on simplices.
Implementation of range helper class.
unsigned int i_eval_point_
index of point in EvalPoint object
unsigned int i_element_
mesh_idx of ElementAccessor appropriate to element
unsigned int i_reg_
region_idx of element
Helper struct stores data for initizalize descentants of FieldAlgorithmBase.
#define END_TIMER(tag)
Ends a timer with specified tag.
#define ADD_CALLS(n_calls)
Increase number of calls in actual timer.
#define START_TIMER(tag)
Starts a timer with specified tag.
@ update_values
Shape function values.