v0.15.0
Loading...
Searching...
No Matches
NormsOperators.cpp
Go to the documentation of this file.
1/** \file NormsOperators.cpp
2
3\brief User data operators for calculating norms and differences between
4fields
5
6*/
7
8namespace MoFEM {
9
11 boost::shared_ptr<VectorDouble> data_ptr, SmartPetscObj<Vec> data_vec,
12 const int index, boost::shared_ptr<VectorDouble> diff_data_ptr,
13 boost::shared_ptr<Range> ent_ptr)
16 dataPtr(data_ptr), dataVec(data_vec), iNdex(index),
17 diffDataPtr(diff_data_ptr), entsPtr(ent_ptr) {
18 if (!dataPtr)
19 THROW_MESSAGE("Pointer is not set");
20 if (!diffDataPtr)
22}
23
27
28 // check if entity is in the range
29 if (entsPtr) {
30 if (entsPtr->find(this->getFEEntityHandle()) == entsPtr->end())
32 }
33
34 // calculate the difference between data pointers and save them in diffDataPtr
35 if (dataPtr != diffDataPtr)
37
38 // get number of integration points
39 const auto nb_integration_points = getGaussPts().size2();
40 // get element volume
41 const double vol = getMeasure();
42 // get integration weights
44 // get values
45 auto t_data = getFTensor0FromVec(*diffDataPtr);
46 // initialise double to store norm values
47 double norm_on_element = 0.;
48 // loop over integration points
49 for (int gg = 0; gg != nb_integration_points; gg++) {
50 // add to element norm
51 norm_on_element += t_w * t_data * t_data;
52 // move to another integration weight
53 ++t_w;
54 // move to another data values
55 ++t_data;
56 }
57 // scale with volume of the element
58 norm_on_element *= vol;
59 // add to dataVec at iNdex position
60 CHKERR VecSetValue(dataVec, iNdex, norm_on_element, ADD_VALUES);
61
63}
64
65template <int DIM>
67 boost::shared_ptr<MatrixDouble> data_ptr, SmartPetscObj<Vec> data_vec,
68 const int index, boost::shared_ptr<MatrixDouble> diff_data_ptr,
69 boost::shared_ptr<Range> ent_ptr)
72 dataPtr(data_ptr), dataVec(data_vec), iNdex(index),
73 diffDataPtr(diff_data_ptr), entsPtr(ent_ptr) {
74 if (!dataPtr)
75 THROW_MESSAGE("Pointer is not set");
76 if (!diffDataPtr)
78}
79
80template <int DIM>
82OpCalcNormL2Tensor1<DIM>::doWork(int side, EntityType type,
85
86 // check if entity is in the range
87 if (entsPtr) {
88 if (entsPtr->find(this->getFEEntityHandle()) == entsPtr->end())
90 }
91
92 // calculate the difference between data pointers and save them in diffDataPtr
93 if (dataPtr != diffDataPtr)
94 *diffDataPtr -= *dataPtr;
95
96 // Declare FTensor index
97 FTensor::Index<'i', DIM> i;
98 // get number of integration points
99 const auto nb_integration_points = getGaussPts().size2();
100 // get element volume
101 const double vol = getMeasure();
102 // get integration weights
103 auto t_w = getFTensor0IntegrationWeight();
104 // get vector values
105 auto t_data = getFTensor1FromMat<DIM>(*diffDataPtr);
106 // initialise double to store norm values
107 double norm_on_element = 0.;
108 // loop over integration points
109 for (int gg = 0; gg != nb_integration_points; ++gg) {
110 // add to element norm
111 norm_on_element += t_w * (t_data(i) * t_data(i));
112 // move to another integration weight
113 ++t_w;
114 // move to another data values
115 ++t_data;
116 }
117 // scale with volume of the element
118 norm_on_element *= vol;
119 // add to dataVec at iNdex position
120 CHKERR VecSetValue(dataVec, iNdex, norm_on_element, ADD_VALUES);
121
123}
124
125template <int DIM_1, int DIM_2>
127 boost::shared_ptr<MatrixDouble> data_ptr, SmartPetscObj<Vec> data_vec,
128 const int index, boost::shared_ptr<MatrixDouble> diff_data_ptr,
129 boost::shared_ptr<Range> ent_ptr)
132 dataPtr(data_ptr), dataVec(data_vec), iNdex(index),
133 diffDataPtr(diff_data_ptr), entsPtr(ent_ptr) {
134 if (!dataPtr)
135 THROW_MESSAGE("Pointer is not set");
136 if (!diffDataPtr)
138}
139
140template <int DIM_1, int DIM_2>
145
146 // check if entity is in the range
147 if (entsPtr) {
148 if (entsPtr->find(this->getFEEntityHandle()) == entsPtr->end())
150 }
151
152 // calculate the difference between data pointers and save them in diffDataPtr
153 if (dataPtr != diffDataPtr)
154 *diffDataPtr -= *dataPtr;
155
156 // Declare FTensor index
157 FTensor::Index<'i', DIM_1> i;
158 FTensor::Index<'j', DIM_2> j;
159 // get number of integration points
160 const auto nb_integration_points = getGaussPts().size2();
161 // get element volume
162 const double vol = getMeasure();
163 // get integration weights
164 auto t_w = getFTensor0IntegrationWeight();
165 // get vector values
166 auto t_data = getFTensor2FromMat<DIM_1, DIM_2>(*diffDataPtr);
167 // initialise double to store norm values
168 double norm_on_element = 0.;
169 // loop over integration points
170 for (int gg = 0; gg != nb_integration_points; gg++) {
171 // add to element norm
172 norm_on_element += t_w * (t_data(i, j) * t_data(i, j));
173 // move to another integration weight
174 ++t_w;
175 // move to another data values
176 ++t_data;
177 }
178 // scale with volume of the element
179 norm_on_element *= vol;
180 // add to dataVec at iNdex position
181 CHKERR VecSetValue(dataVec, iNdex, norm_on_element, ADD_VALUES);
182
184}
185
187 boost::shared_ptr<VectorDouble> data_ptr, ScalarFun scalar_function)
190 dataPtr(data_ptr), sFunc(scalar_function) {}
191
195
196 // get number of integration points
197 const int nb_integration_pts = getGaussPts().size2();
198 // resize dataPtr to store values at integration points
199 dataPtr->resize(nb_integration_pts);
200 // get values
201 auto t_val = getFTensor0FromVec(*dataPtr);
202
203 // get coordinates at integration point
204 auto t_coords = getFTensor1CoordsAtGaussPts();
205 // loop over integration points
206 for (int gg = 0; gg != nb_integration_pts; ++gg) {
207 // set value at integration point to value of scalar function at the
208 // coordinates
209 t_val = sFunc(t_coords(0), t_coords(1), t_coords(2));
210 // move to another value
211 ++t_val;
212 // move to another coordinates
213 ++t_coords;
214 }
215
217}
218
219template <int SPACE_DIM, int BASE_DIM>
221 boost::shared_ptr<MatrixDouble> data_ptr, VectorFunc vector_function)
224 dataPtr(data_ptr), vFunc(vector_function) {}
225
226template <int SPACE_DIM, int BASE_DIM>
228 int side, EntityType type, EntitiesFieldData::EntData &data) {
230
231 // Declare FTensor index set to space dimension
233 // get number of integration points
234 const int nb_integration_pts = getGaussPts().size2();
235 // resize dataPtr to store values at integration points
236 dataPtr->resize(BASE_DIM, nb_integration_pts);
237 // get values
238 auto t_val = getFTensor1FromMat<BASE_DIM>(*(dataPtr));
239
240 // get coordinates at integration point
241 auto t_coords = getFTensor1CoordsAtGaussPts();
242 // loop over integration points
243 for (int gg = 0; gg != nb_integration_pts; ++gg) {
244 // get function values of vector function at the coordinates
245 auto func_val = vFunc(t_coords(0), t_coords(1), t_coords(2));
246 // translate function values to tensor
247 auto t_func_val = getFTensor1FromArray<SPACE_DIM, SPACE_DIM>(func_val);
248 // set values at integration point to function values at the coordinates
249 t_val(i) = t_func_val(i);
250
251 // move to another value
252 ++t_val;
253 // move to another coordinates
254 ++t_coords;
255 }
256
258}
259
260template struct OpCalcNormL2Tensor1<2>;
261template struct OpCalcNormL2Tensor1<3>;
262template struct OpCalcNormL2Tensor2<2, 2>;
263template struct OpCalcNormL2Tensor2<3, 3>;
264template struct OpGetTensor1fromFunc<2, 2>;
265template struct OpGetTensor1fromFunc<2, 3>;
266template struct OpGetTensor1fromFunc<3, 3>;
267
268} // namespace MoFEM
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr int SPACE_DIM
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ NOSPACE
Definition definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
constexpr int BASE_DIM
boost::function< double(const double, const double, const double)> ScalarFun
Scalar function type.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
boost::function< VectorDouble(const double, const double, const double)> VectorFunc
auto getFTensor1FromArray(VectorDouble &data)
Get FTensor1 from array.
Data on single entity (This is passed as argument to DataOperator::doWork)
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
structure to get information form mofem into EntitiesFieldData
boost::shared_ptr< VectorDouble > dataPtr
boost::shared_ptr< VectorDouble > diffDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate norm of scalar values at integration points
OpCalcNormL2Tensor0(boost::shared_ptr< VectorDouble > data_ptr, SmartPetscObj< Vec > data_vec, const int index, boost::shared_ptr< VectorDouble > diff_data_ptr=nullptr, boost::shared_ptr< Range > ent_ptr=nullptr)
SmartPetscObj< Vec > dataVec
boost::shared_ptr< Range > entsPtr
Get norm of input MatrixDouble for Tensor1.
boost::shared_ptr< MatrixDouble > diffDataPtr
OpCalcNormL2Tensor1(boost::shared_ptr< MatrixDouble > data_ptr, SmartPetscObj< Vec > data_vec, const int index, boost::shared_ptr< MatrixDouble > diff_data_ptr=nullptr, boost::shared_ptr< Range > ent_ptr=nullptr)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate norm of vector values at integration points
boost::shared_ptr< MatrixDouble > dataPtr
Get norm of input MatrixDouble for Tensor2.
OpCalcNormL2Tensor2(boost::shared_ptr< MatrixDouble > data_ptr, SmartPetscObj< Vec > data_vec, const int index, boost::shared_ptr< MatrixDouble > diff_data_ptr=nullptr, boost::shared_ptr< Range > ent_ptr=nullptr)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate norm of tensor values at integration points
boost::shared_ptr< MatrixDouble > diffDataPtr
boost::shared_ptr< MatrixDouble > dataPtr
OpGetTensor0fromFunc(boost::shared_ptr< VectorDouble > data_ptr, ScalarFun scalar_function)
boost::shared_ptr< VectorDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of scalar function at integration points
OpGetTensor1fromFunc(boost::shared_ptr< MatrixDouble > data_ptr, VectorFunc vector_function)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of vector function at integration points
intrusive_ptr for managing petsc objects