v0.14.0
DispMap.hpp
Go to the documentation of this file.
1 /* This file is part of MoFEM.
2 
3  * MoFEM is free software: you can redistribute it and/or modify it under
4  * the terms of the GNU Lesser General Public License as published by the
5  * Free Software Foundation, either version 3 of the License, or (at your
6  * option) any later version.
7  *
8  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
9  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
11  * License for more details.
12  *
13  * You should have received a copy of the GNU Lesser General Public
14  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
15 
16 #ifndef __DISP_MAP_HPP__
17 #define __DISP_MAP_HPP__
18 
19 namespace CellEngineering {
20 
21 struct DispMapFe : public MoFEM::FaceElementForcesAndSourcesCore {
22  DispMapFe(MoFEM::Interface &mField)
23  : MoFEM::FaceElementForcesAndSourcesCore(mField) {}
24  /** \brief Set integrate rule
25  */
26  int getRule(int order) { return 2 * order + 1; };
27 };
28 
29 struct DataFromFiles {
30 
31  string &fIle;
32 
33  double sCale;
34  double cUbe_size;
35  vector<vector<double>> main_vector;
36  double lAmbda;
37 
38  DataFromFiles(string &file) : fIle(file) {}
39 
40  MoFEMErrorCode loadFileData() {
41  MoFEMFunctionBegin;
42  ifstream in(fIle.c_str());
43  if (!in.is_open()) {
44  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
45  "file with data points not found");
46  // cout << ">> couldn't read the file <<"<< endl;
47  }
48 
49  typedef boost::tokenizer<boost::escaped_list_separator<char>> Tokenizer;
50 
51  string line;
52 
53  main_vector.clear();
54  while (getline(in, line)) {
55  vector<double> vec;
56  Tokenizer tok(line);
57  for (Tokenizer::iterator it(tok.begin()), end(tok.end()); it != end;
58  ++it) {
59  vec.push_back(atof(it->c_str()));
60  }
61 
62  main_vector.push_back(vec);
63  }
64  reverse(main_vector.begin(), main_vector.end());
65  MoFEMFunctionReturn(0);
66  }
67 
68  MoFEMErrorCode fromOptions() {
69  MoFEMFunctionBegin;
70  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Get parameters", "none");
71 
72  cUbe_size = 4;
73  CHKERR PetscOptionsScalar("-cube_size", "get cube size for the averaging",
74  "", 0, &cUbe_size, PETSC_NULL);
75 
76  sCale = 1;
77  CHKERR PetscOptionsScalar(
78  "-scale", "scale the distance between voxels (eg. from mm to m)", "",
79  sCale, &sCale, PETSC_NULL);
80 
81  ierr = PetscOptionsEnd();
82  CHKERRQ(ierr);
83 
84  MoFEMFunctionReturn(0);
85  }
86 
87  /**
88  * returns vertors of data for given index of the point
89  * @param vec_idx reference vector of indices
90  * @param vec_data reference vector of data (colors)
91  * @return error code
92  */
93 
94  double getDataForGivenCoordinate(const double x, const double y,
95  vector<vector<double>> &main_vector,
96  const double sCale, const double cUbe_size) {
97 
98  const int size_y = main_vector.size(); // flip
99  const int size_x = main_vector.begin()->size();
100  double data = 0;
101 
102  vector<int> vec_ix;
103  vector<int> vec_iy;
104  vector<double> vec_data;
105  vector<double> vec_dist;
106 
107  int nx = ceil(cUbe_size);
108  int ny = ceil(cUbe_size);
109 
110  int ix = ceil(x / sCale);
111  int iy = ceil(y / sCale);
112 
113  vec_ix.resize(2 * nx);
114  vec_iy.resize(2 * ny);
115 
116  int ii = 0;
117  for (int i = 0; i < 2 * nx; i++) {
118  int idx = ix - nx + i;
119  if (idx >= size_x / 2 || idx < -size_x / 2)
120  continue;
121  vec_ix[ii++] = idx;
122  }
123  vec_ix.resize(ii);
124  ii = 0;
125  for (int i = 0; i < 2 * ny; i++) {
126  int idx = iy - ny + i;
127  if (idx >= size_y / 2 || idx < -size_y / 2)
128  continue;
129  vec_iy[ii++] = idx;
130  }
131  vec_iy.resize(ii);
132  ii = 0;
133 
134  vec_data.resize(vec_iy.size() * vec_ix.size());
135  vec_dist.resize(vec_iy.size() * vec_ix.size());
136 
137  for (vector<int>::iterator it_iy = vec_iy.begin(); it_iy != vec_iy.end();
138  it_iy++) {
139  for (vector<int>::iterator it_ix = vec_ix.begin(); it_ix != vec_ix.end();
140  it_ix++) {
141  vec_data[ii] =
142  main_vector[*it_iy + size_y / 2][*it_ix + size_x / 2]; // flip
143  vec_dist[ii] = ((*it_iy * sCale - y) * (*it_iy * sCale - y) +
144  (*it_ix * sCale - x) * (*it_ix * sCale - x));
145 
146  ii++;
147  }
148  }
149 
150  // gaussian smoothing
151  vector<double> kernel;
152  kernel.resize(vec_data.size());
153  int i = 0;
154  double sigma = 10;
155  double sum = 0;
156  const double m = (1 / sqrt(M_PI * 2 * sigma * sigma));
157  for (int ii = 0; ii < vec_dist.size(); ii++) {
158  kernel[i] = m * exp(-(vec_dist[ii]) /
159  (2 * sigma * sigma)); // distance is already squared
160  sum += kernel[i++];
161  }
162  ii = 0;
163  for (vector<double>::iterator vec_itr = vec_data.begin();
164  vec_itr != vec_data.end(); vec_itr++) {
165  kernel[ii] /= sum;
166  data += (*vec_itr) * kernel[ii++];
167  }
168  return data;
169  }
170 
171 };
172 
173 struct CommonData {
174 
175  VectorDouble dispAtGaussPts; ///< Values at integration Pts
176  VectorDouble locF; ///< Local element rhs vector
177  MatrixDouble locA; ///< Local element matrix
178  MatrixDouble transLocA;
179  Mat globA; ///< Global matrix
180  Vec globF;
181  MatrixDouble invJac;
182  double lAmbda;
183 
184  MoFEMErrorCode getParameters() {
185 
186  MoFEMFunctionBegin;
187  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "none");
188 
189  this->lAmbda = 1;
190  CHKERR PetscOptionsScalar("-lambda", "lambda parameter", "", 1, &lAmbda,
191  PETSC_NULL);
192  ierr = PetscOptionsEnd();
193  CHKERRG(ierr);
194 
195  MoFEMFunctionReturn(0);
196  }
197 };
198 
199 struct OpCalculateLhs
200  : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
201 
202  CommonData &commonData;
203  DataFromFiles &dataFromFiles;
204 
205  OpCalculateLhs(const std::string &row_field, const std::string &col_field,
206  CommonData &common_data, DataFromFiles &data_from_files)
207  : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
208  row_field, col_field, UserDataOperator::OPROWCOL),
209  commonData(common_data), dataFromFiles(data_from_files) {
210  sYmm = true;
211  }
212 
213  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
214  EntityType col_type,
215  DataForcesAndSourcesCore::EntData &row_data,
216  DataForcesAndSourcesCore::EntData &col_data) {
217 
218  MoFEMFunctionBegin;
219 
220  int nb_rows = row_data.getIndices().size(); // number of rows
221  int nb_cols = col_data.getIndices().size(); // number of columns
222  if (nb_rows == 0)
223  MoFEMFunctionReturnHot(0);
224  if (nb_cols == 0)
225  MoFEMFunctionReturnHot(0);
226 
227  commonData.locA.resize(nb_rows, nb_cols, false);
228  commonData.locA.clear();
229  commonData.getParameters();
230  int nb_gauss_pts = row_data.getN().size1();
231  for (int gg = 0; gg < nb_gauss_pts; gg++) {
232  double area = getArea() * getGaussPts()(2, gg);
233  // if(getHoGaussPtsDetJac().size()>0) {
234  // area *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
235  // }
236  VectorAdaptor row_base_function = row_data.getN(gg);
237  VectorAdaptor col_base_function = col_data.getN(gg);
238  commonData.locA +=
239  area * outer_prod(row_base_function, col_base_function);
240  MatrixAdaptor row_diff_base_function = row_data.getDiffN(gg);
241  MatrixAdaptor col_diff_base_function = col_data.getDiffN(gg);
242 
243  // const double lambda = 1;
244 
245  commonData.locA +=
246  area * commonData.lAmbda * dataFromFiles.sCale * dataFromFiles.sCale *
247  prod(row_diff_base_function, trans(col_diff_base_function));
248  }
249  CHKERR MatSetValues(commonData.globA, row_data.getIndices().size(),
250  &row_data.getIndices()[0], col_data.getIndices().size(),
251  &col_data.getIndices()[0], &commonData.locA(0, 0),
252  ADD_VALUES);
253 
254  // Assemble off symmetric side
255  if (row_side != col_side || row_type != col_type) {
256  commonData.transLocA.resize(commonData.locA.size2(),
257  commonData.locA.size1(), false);
258  noalias(commonData.transLocA) = trans(commonData.locA);
259  CHKERR MatSetValues(
260  commonData.globA, col_data.getIndices().size(),
261  &col_data.getIndices()[0], row_data.getIndices().size(),
262  &row_data.getIndices()[0], &commonData.transLocA(0, 0), ADD_VALUES);
263  }
264 
265  MoFEMFunctionReturn(0);
266  }
267 };
268 
269 struct OpCalulatefRhoAtGaussPts
270  : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
271 
272  CommonData &commonData;
273  DataFromFiles &dataFromFiles;
274 
275  OpCalulatefRhoAtGaussPts(const string &row_field, CommonData &common_data,
276  DataFromFiles &data_from_files)
277  : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
278  row_field, row_field, UserDataOperator::OPROW),
279  commonData(common_data), dataFromFiles(data_from_files) {}
280 
281  MoFEMErrorCode doWork(int row_side, EntityType row_type,
282  DataForcesAndSourcesCore::EntData &row_data) {
283  MoFEMFunctionBegin;
284 
285  if (row_type != MBVERTEX)
286  MoFEMFunctionReturnHot(0);
287 
288  int nb_rows = row_data.getIndices().size();
289  if (nb_rows == 0)
290  MoFEMFunctionReturnHot(0);
291 
292  FTensor::Index<'i', 3> i;
293 
294  FTensor::Tensor1<double *, 3> t1_xyz(&getCoordsAtGaussPts()(0, 0),
295  &getCoordsAtGaussPts()(0, 1),
296  &getCoordsAtGaussPts()(0, 2), 3);
297  const int nb_gauss_pts = row_data.getN().size1();
298  commonData.dispAtGaussPts.resize(nb_gauss_pts, false);
299  for (int gg = 0; gg < nb_gauss_pts; gg++) {
300 
301  commonData.dispAtGaussPts[gg] = dataFromFiles.getDataForGivenCoordinate(
302  t1_xyz(0), t1_xyz(1), dataFromFiles.main_vector, dataFromFiles.sCale,
303  dataFromFiles.cUbe_size);
304  ++t1_xyz;
305  }
306 
307  MoFEMFunctionReturn(0);
308  }
309 };
310 
311 struct OpAssmbleRhs
312  : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
313 
314  CommonData &commonData;
315 
316  OpAssmbleRhs(const string &row_field, CommonData &common_data)
317  : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
318  row_field, row_field, UserDataOperator::OPROW),
319  commonData(common_data) {}
320 
321  MoFEMErrorCode doWork(int row_side, EntityType row_type,
322  DataForcesAndSourcesCore::EntData &row_data) {
323  MoFEMFunctionBegin;
324 
325  int nb_rows = row_data.getIndices().size();
326  if (nb_rows == 0)
327  MoFEMFunctionReturnHot(0);
328  commonData.locF.resize(nb_rows, false);
329  commonData.locF.clear();
330 
331  int nb_gauss_pts = row_data.getN().size1();
332  for (int gg = 0; gg < nb_gauss_pts; gg++) {
333 
334  double area = getArea() * getGaussPts()(2, gg);
335  VectorAdaptor base_function = row_data.getN(gg);
336  commonData.locF += area * commonData.dispAtGaussPts[gg] * base_function;
337  }
338 
339  CHKERR VecSetValues(commonData.globF, row_data.getIndices().size(),
340  &row_data.getIndices()[0], &commonData.locF[0],
341  ADD_VALUES);
342 
343  MoFEMFunctionReturn(0);
344  }
345 };
346 
347 } // namespace CellEngineering
348 
349 #endif //__DISP_MAP_HPP__
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
CellEngineering::DataFromFiles
Definition: DispMap.hpp:29
CellEngineering
Definition: cell_forces.cpp:412
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1644
CellEngineering::DataFromFiles::cUbe_size
double cUbe_size
Definition: DispMap.hpp:34
CellEngineering::DataFromFiles::DataFromFiles
DataFromFiles(string &file)
Definition: DispMap.hpp:38
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
CellEngineering::OpCalculateLhs::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: DispMap.hpp:213
CellEngineering::CommonData::globF
Vec globF
Definition: DispMap.hpp:180
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
CellEngineering::CommonData::locA
MatrixDouble locA
Local element matrix.
Definition: DispMap.hpp:177
CellEngineering::DataFromFiles::fromOptions
MoFEMErrorCode fromOptions()
Definition: DispMap.hpp:68
CellEngineering::OpCalculateLhs::commonData
CommonData & commonData
Definition: DispMap.hpp:202
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROWCOL
@ OPROWCOL
operator doWork is executed on FE rows &columns
Definition: ForcesAndSourcesCore.hpp:569
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::ForcesAndSourcesCore::UserDataOperator::getGaussPts
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Definition: ForcesAndSourcesCore.hpp:1236
CellEngineering::CommonData::invJac
MatrixDouble invJac
Definition: DispMap.hpp:181
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
CellEngineering::DispMapFe
Definition: DispMap.hpp:21
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
CellEngineering::OpCalculateLhs::OpCalculateLhs
OpCalculateLhs(const std::string &row_field, const std::string &col_field, CommonData &common_data, DataFromFiles &data_from_files)
Definition: DispMap.hpp:205
CellEngineering::OpCalculateLhs
Definition: DispMap.hpp:199
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
CellEngineering::DataFromFiles::main_vector
vector< vector< double > > main_vector
Definition: DispMap.hpp:35
CellEngineering::OpAssmbleRhs::OpAssmbleRhs
OpAssmbleRhs(const string &row_field, CommonData &common_data)
Definition: DispMap.hpp:316
CellEngineering::CommonData::lAmbda
double lAmbda
Definition: DispMap.hpp:182
CellEngineering::OpAssmbleRhs::commonData
CommonData & commonData
Definition: DispMap.hpp:314
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
CellEngineering::CommonData::globA
Mat globA
Global matrix.
Definition: DispMap.hpp:179
CellEngineering::OpCalulatefRhoAtGaussPts::OpCalulatefRhoAtGaussPts
OpCalulatefRhoAtGaussPts(const string &row_field, CommonData &common_data, DataFromFiles &data_from_files)
Definition: DispMap.hpp:275
MoFEM::ForcesAndSourcesCore::UserDataOperator::getCoordsAtGaussPts
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
Definition: ForcesAndSourcesCore.hpp:1265
MoFEM::Types::VectorAdaptor
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:115
CellEngineering::DataFromFiles::fIle
string & fIle
Definition: DispMap.hpp:31
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
CellEngineering::DispMapFe::getRule
int getRule(int order)
Set integrate rule.
Definition: DispMap.hpp:26
FaceElementForcesAndSourcesCore
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
FTensor::Index< 'i', 3 >
CellEngineering::OpCalulatefRhoAtGaussPts
Definition: DispMap.hpp:269
CellEngineering::OpCalculateLhs::dataFromFiles
DataFromFiles & dataFromFiles
Definition: DispMap.hpp:203
CellEngineering::CommonData::locF
VectorDouble locF
Local element rhs vector.
Definition: DispMap.hpp:176
CellEngineering::DataFromFiles::sCale
double sCale
Definition: DispMap.hpp:33
CellEngineering::OpCalulatefRhoAtGaussPts::commonData
CommonData & commonData
Definition: DispMap.hpp:272
CellEngineering::OpCalulatefRhoAtGaussPts::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: DispMap.hpp:281
CellEngineering::CommonData::getParameters
MoFEMErrorCode getParameters()
Definition: DispMap.hpp:184
CommonData
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:22
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
CellEngineering::DispMapFe::DispMapFe
DispMapFe(MoFEM::Interface &mField)
Definition: DispMap.hpp:22
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
CellEngineering::DataFromFiles::loadFileData
MoFEMErrorCode loadFileData()
Definition: DispMap.hpp:40
CellEngineering::CommonData
Definition: CellForces.hpp:31
CellEngineering::OpAssmbleRhs
Definition: DispMap.hpp:311
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
CellEngineering::CommonData::dispAtGaussPts
VectorDouble dispAtGaussPts
Values at integration Pts.
Definition: DispMap.hpp:175
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::Types::MatrixAdaptor
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
Definition: Types.hpp:132
CellEngineering::CommonData::transLocA
MatrixDouble transLocA
Definition: DispMap.hpp:178
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator::getArea
double getArea()
get area of face
Definition: FaceElementForcesAndSourcesCore.hpp:239
CellEngineering::OpAssmbleRhs::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: DispMap.hpp:321
CellEngineering::DataFromFiles::lAmbda
double lAmbda
Definition: DispMap.hpp:36
CellEngineering::OpCalulatefRhoAtGaussPts::dataFromFiles
DataFromFiles & dataFromFiles
Definition: DispMap.hpp:273
MoFEM::DataOperator::sYmm
bool sYmm
If true assume that matrix is symmetric structure.
Definition: DataOperators.hpp:82
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
CellEngineering::DataFromFiles::getDataForGivenCoordinate
double getDataForGivenCoordinate(const double x, const double y, vector< vector< double >> &main_vector, const double sCale, const double cUbe_size)
Definition: DispMap.hpp:94
MoFEM::ForcesAndSourcesCore::mField
Interface & mField
Definition: ForcesAndSourcesCore.hpp:24
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
operator doWork function is executed on FE rows
Definition: ForcesAndSourcesCore.hpp:567
Generated by Doxygen 1.8.17 and hosted at University of Glasgow