LifeV
RBFrescaledScalar.hpp
Go to the documentation of this file.
1 //@HEADER
2 /*
3 *******************************************************************************
4 
5  Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
6  Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
7 
8  This file is part of LifeV.
9 
10  LifeV is free software; you can redistribute it and/or modify
11  it under the terms of the GNU Lesser General Public License as published by
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  LifeV is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public License
21  along with LifeV. If not, see <http://www.gnu.org/licenses/>.
22 
23 *******************************************************************************
24 */
25 //@HEADER
26 
27 /*!
28  @file
29  @brief A short description of the file content
30 
31  @author Davide Forti <forti@mathicsepc48.epfl.ch>
32  @date 14 Mar 2013
33 
34  A more detailed description of the file (if necessary)
35  */
36 
37 #ifndef RBFRESCALEDSCALAR_H
38 #define RBFRESCALEDSCALAR_H 1
39 
40 #include <lifev/core/interpolation/RBFInterpolation.hpp>
41 
42 namespace LifeV {
43 
44 template <typename mesh_Type>
45 class RBFrescaledScalar: public RBFInterpolation<mesh_Type>
46 {
47 public:
48 
50 
51  typedef VectorEpetra vector_Type;
53 
54  typedef MatrixEpetra<double> matrix_Type;
56 
57  typedef std::vector<int> flagContainer_Type;
58 
60 
63 
64  typedef GhostHandler<mesh_Type> neighbors_Type;
66 
67  typedef LifeV::Preconditioner basePrec_Type;
69 
70  typedef LifeV::PreconditionerIfpack prec_Type;
72 
74 
76 
77  virtual ~RBFrescaledScalar();
78 
79  void setup ( meshPtr_Type fullMeshKnown, meshPtr_Type localMeshKnown, meshPtr_Type fullMeshUnknown, meshPtr_Type localMeshUnknown, flagContainer_Type flags );
80 
81  void setupRBFData (vectorPtr_Type KnownField, vectorPtr_Type UnknownField, GetPot datafile, parameterList_Type belosList);
82 
83  void buildOperators();
84 
85  void interpolationOperator();
86 
87  void projectionOperator();
88 
89  void buildRhs();
90 
92 
93  void identifyNodes (meshPtr_Type LocalMesh, std::unordered_set<ID>& GID_nodes, vectorPtr_Type CheckVector);
94 
95  bool isInside (ID pointMarker, flagContainer_Type Flags);
96 
97  double rbf (double x1, double y1, double z1, double x2, double y2, double z2, double radius);
98 
99  void interpolate();
100 
101  void solution (vectorPtr_Type& Solution);
102 
103  void solutionrbf (vectorPtr_Type& Solution_rbf);
104 
105  void updateRhs(vectorPtr_Type newRhs);
106 
107  void setRadius ( double radius );
108 
109 private:
110 
131  double M_radius;
132 };
133 
134 template <typename mesh_Type>
136  M_radius(0)
137 {}
138 
139 template <typename mesh_Type>
141 {}
142 
143 template <typename mesh_Type>
144 void RBFrescaledScalar<mesh_Type>::setup ( meshPtr_Type fullMeshKnown, meshPtr_Type localMeshKnown, meshPtr_Type fullMeshUnknown, meshPtr_Type localMeshUnknown, flagContainer_Type flags )
145 {
146  M_fullMeshKnown = fullMeshKnown;
147  M_localMeshKnown = localMeshKnown;
148  M_fullMeshUnknown = fullMeshUnknown;
149  M_localMeshUnknown = localMeshUnknown;
150  M_flags = flags;
151 }
152 
153 template <typename mesh_Type>
154 void RBFrescaledScalar<mesh_Type>::setRadius ( double radius )
155 {
156  M_radius = radius;
157 }
158 
159 template <typename mesh_Type>
160 void RBFrescaledScalar<mesh_Type>::setupRBFData (vectorPtr_Type KnownField, vectorPtr_Type UnknownField, GetPot datafile, parameterList_Type belosList)
161 {
162  M_knownField = KnownField;
163  M_unknownField = UnknownField;
164  M_datafile = datafile;
165  M_belosList = belosList;
166 }
167 
168 
169 template <typename mesh_Type>
171 {
172  LifeChrono TimeBuilding;
173  TimeBuilding.start();
174  this->interpolationOperator();
175  this->projectionOperator();
176  TimeBuilding.stop();
177  if(M_knownField->mapPtr()->commPtr()->MyPID()==0)
178  std::cout << "Time to assembly operators = " << TimeBuilding.diff() << std::endl;
179 
180  this->buildRhs();
181  this->interpolateCostantField();
182 }
183 
184 template <typename mesh_Type>
186 {
187  ASSERT(M_radius!=0, "Please set the basis radius using RBFrescaledScalar<mesh_Type>::setRadius(double radius)")
188  this->identifyNodes (M_localMeshKnown, M_GIdsKnownMesh, M_knownField);
189  M_neighbors.reset ( new neighbors_Type ( M_fullMeshKnown, M_localMeshKnown, M_knownField->mapPtr(), M_knownField->mapPtr()->commPtr() ) );
190  if (M_flags[0] == -1)
191  {
192  M_neighbors->setUpNeighbors();
193  }
194  else
195  {
196  M_neighbors->createPointPointNeighborsList (M_flags);
197  }
198 
199  int LocalNodesNumber = M_GIdsKnownMesh.size();
200 
201  std::vector<std::unordered_set<ID> > MatrixGraph (LocalNodesNumber);
202  int* ElementsPerRow = new int[LocalNodesNumber];
203  int* GlobalID = new int[LocalNodesNumber];
204  int k = 0;
205  int Max_entries = 0;
206 
207  for (std::unordered_set<ID>::iterator it = M_GIdsKnownMesh.begin(); it != M_GIdsKnownMesh.end(); ++it)
208  {
209  GlobalID[k] = *it;
210  MatrixGraph[k] = M_neighbors->neighborsWithinRadius (GlobalID[k], M_radius);
211  MatrixGraph[k].insert (GlobalID[k]);
212  ElementsPerRow[k] = MatrixGraph[k].size();
213  if (ElementsPerRow[k] > Max_entries)
214  {
215  Max_entries = ElementsPerRow[k];
216  }
217  ++k;
218  }
219 
220  M_interpolationOperatorMap.reset (new map_Type (-1, LocalNodesNumber, GlobalID, M_knownField->mapPtr()->commPtr() ) );
221  M_interpolationOperator.reset (new matrix_Type (*M_interpolationOperatorMap, ElementsPerRow) );
222 
223  int* Indices = new int[Max_entries];
224  double* Values = new double[Max_entries];
225 
226  for ( int i = 0 ; i < LocalNodesNumber; ++i )
227  {
228  k = 0;
229  for ( std::unordered_set<ID>::iterator it = MatrixGraph[i].begin(); it != MatrixGraph[i].end(); ++it)
230  {
231  Indices[k] = *it;
232  Values[k] = rbf ( M_fullMeshKnown->point (GlobalID[i]).x(),
233  M_fullMeshKnown->point (GlobalID[i]).y(),
234  M_fullMeshKnown->point (GlobalID[i]).z(),
235  M_fullMeshKnown->point (*it).x(),
236  M_fullMeshKnown->point (*it).y(),
237  M_fullMeshKnown->point (*it).z(),
238  M_radius);
239  ++k;
240  }
241  M_interpolationOperator->matrixPtr()->InsertGlobalValues (GlobalID[i], k, Values, Indices);
242  }
243  M_interpolationOperator->globalAssemble();
244  delete[] Indices;
245  delete[] Values;
246  delete[] ElementsPerRow;
247  delete[] GlobalID;
248 }
249 
250 template <typename mesh_Type>
252 {
253 
254  this->identifyNodes (M_localMeshUnknown, M_GIdsUnknownMesh, M_unknownField);
255 
256  int LocalNodesNumber = M_GIdsUnknownMesh.size();
257 
258  std::vector<std::unordered_set<ID> > MatrixGraph (LocalNodesNumber);
259  int* ElementsPerRow = new int[LocalNodesNumber];
260  int* GlobalID = new int[LocalNodesNumber];
261  int k = 0;
262  int Max_entries = 0;
263  double d;
264  double d_min;
265  int nearestPoint;
266 
267  for (std::unordered_set<ID>::iterator it = M_GIdsUnknownMesh.begin(); it != M_GIdsUnknownMesh.end(); ++it)
268  {
269  GlobalID[k] = *it;
270  d_min = 100;
271  for (int j = 0; j < M_fullMeshKnown->numVertices(); ++j)
272  {
273  if ( M_flags[0] == -1 || this->isInside (M_fullMeshKnown->point (j).markerID(), M_flags) )
274  {
275  d = std::sqrt ( pow (M_fullMeshKnown->point (j).x() - M_fullMeshUnknown->point (GlobalID[k]).x(), 2)
276  + pow (M_fullMeshKnown->point (j).y() - M_fullMeshUnknown->point (GlobalID[k]).y(), 2)
277  + pow (M_fullMeshKnown->point (j).z() - M_fullMeshUnknown->point (GlobalID[k]).z(), 2) );
278  if (d < d_min)
279  {
280  d_min = d;
281  nearestPoint = M_fullMeshKnown->point (j).id();
282  }
283  }
284  }
285 
286  MatrixGraph[k] = M_neighbors->neighborsWithinRadius (nearestPoint, M_radius);
287  MatrixGraph[k].insert (nearestPoint);
288  ElementsPerRow[k] = MatrixGraph[k].size();
289  if (ElementsPerRow[k] > Max_entries)
290  {
291  Max_entries = ElementsPerRow[k];
292  }
293  ++k;
294  }
295 
296  M_projectionOperatorMap.reset (new map_Type (-1, LocalNodesNumber, GlobalID, M_unknownField->mapPtr()->commPtr() ) );
297  M_projectionOperator.reset (new matrix_Type (*M_projectionOperatorMap, ElementsPerRow) );
298 
299  int* Indices = new int[Max_entries];
300  double* Values = new double[Max_entries];
301 
302  for ( int i = 0 ; i < LocalNodesNumber; ++i )
303  {
304  k = 0;
305  for ( std::unordered_set<ID>::iterator it = MatrixGraph[i].begin(); it != MatrixGraph[i].end(); ++it)
306  {
307  Indices[k] = *it;
308  Values[k] = rbf ( M_fullMeshUnknown->point (GlobalID[i]).x(),
309  M_fullMeshUnknown->point (GlobalID[i]).y(),
310  M_fullMeshUnknown->point (GlobalID[i]).z(),
311  M_fullMeshKnown->point (*it).x(),
312  M_fullMeshKnown->point (*it).y(),
313  M_fullMeshKnown->point (*it).z(),
314  M_radius);
315  ++k;
316  }
317  M_projectionOperator->matrixPtr()->InsertGlobalValues (GlobalID[i], k, Values, Indices);
318  }
319  M_projectionOperator->globalAssemble (M_interpolationOperatorMap, M_projectionOperatorMap);
320  delete[] Indices;
321  delete[] Values;
322  delete[] ElementsPerRow;
323  delete[] GlobalID;
324 }
325 
326 template <typename mesh_Type>
327 void RBFrescaledScalar<mesh_Type>::buildRhs()
328 {
329  M_RhsF.reset (new vector_Type (*M_interpolationOperatorMap) );
330  M_RhsOne.reset (new vector_Type (*M_interpolationOperatorMap) );
331 
332  M_RhsF->subset (*M_knownField, *M_interpolationOperatorMap, 0, 0);
333  *M_RhsOne += 1;
334 }
335 
336 template <typename mesh_Type>
338 {
339  vectorPtr_Type gamma_one;
340  gamma_one.reset (new vector_Type (*M_interpolationOperatorMap) );
341 
342  // Preconditioner
343  prec_Type* precRawPtr;
344  basePrecPtr_Type precPtr;
345  precRawPtr = new prec_Type;
346  precRawPtr->setDataFromGetPot ( M_datafile, "prec" );
347  precPtr.reset ( precRawPtr );
348 
349  LinearSolver solverOne;
350  solverOne.setCommunicator ( M_knownField->mapPtr()->commPtr() );
351  solverOne.setParameters ( *M_belosList );
352  solverOne.setPreconditioner ( precPtr );
353 
354  solverOne.setOperator (M_interpolationOperator);
355  solverOne.setRightHandSide ( M_RhsOne );
356  solverOne.solve ( gamma_one );
357 
358  M_rbf_one.reset (new vector_Type (*M_projectionOperatorMap) );
359  M_projectionOperator->multiply (false, *gamma_one, *M_rbf_one);
360 }
361 
362 template <typename mesh_Type>
363 void RBFrescaledScalar<mesh_Type>::interpolate()
364 {
365  vectorPtr_Type gamma_f;
366  gamma_f.reset (new vector_Type (*M_interpolationOperatorMap) );
367 
368  // Preconditioner
369  prec_Type* precRawPtr;
370  basePrecPtr_Type precPtr;
371  precRawPtr = new prec_Type;
372  precRawPtr->setDataFromGetPot ( M_datafile, "prec" );
373  precPtr.reset ( precRawPtr );
374 
375  LinearSolver solverRBF;
376  solverRBF.setCommunicator ( M_knownField->mapPtr()->commPtr() );
377  solverRBF.setParameters ( *M_belosList );
378  solverRBF.setPreconditioner ( precPtr );
379 
380  solverRBF.setOperator (M_interpolationOperator);
381  solverRBF.setRightHandSide (M_RhsF);
382  solverRBF.solve (gamma_f);
383 
384  vectorPtr_Type rbf_f;
385  rbf_f.reset (new vector_Type (*M_projectionOperatorMap) );
386 
387  vectorPtr_Type solution;
388  solution.reset (new vector_Type (*M_projectionOperatorMap) );
389 
390  M_projectionOperator->multiply (false, *gamma_f, *rbf_f);
391 
392  *solution = *rbf_f;
393  *solution /= *M_rbf_one;
394 
395  M_unknownField_rbf.reset (new vector_Type (M_unknownField->map() ) );
396  M_unknownField_rbf->subset (*rbf_f, *M_projectionOperatorMap, 0, 0);
397 
398  M_unknownField->subset (*solution, *M_projectionOperatorMap, 0, 0);
399 }
400 
401 
402 template <typename mesh_Type>
403 void RBFrescaledScalar<mesh_Type>::identifyNodes (meshPtr_Type LocalMesh, std::unordered_set<ID>& GID_nodes, vectorPtr_Type CheckVector)
404 {
405 
406  if (M_flags[0] == -1)
407  {
408  for ( UInt i = 0; i < LocalMesh->numVertices(); ++i )
409  if (CheckVector->blockMap().LID ( static_cast<EpetraInt_Type> (LocalMesh->point (i).id()) ) != -1)
410  {
411  GID_nodes.insert (LocalMesh->point (i).id() );
412  }
413  }
414  else
415  {
416  for ( UInt i = 0; i < LocalMesh->numVertices(); ++i )
417  if ( this->isInside (LocalMesh->point (i).markerID(), M_flags) )
418  if (CheckVector->blockMap().LID ( static_cast<EpetraInt_Type> (LocalMesh->point (i).id()) ) != -1)
419  {
420  GID_nodes.insert (LocalMesh->point (i).id() );
421  }
422  }
423 
424 }
425 
426 template <typename mesh_Type>
427 bool RBFrescaledScalar<mesh_Type>::isInside (ID pointMarker, flagContainer_Type flags)
428 {
429  int check = 0;
430  for (UInt i = 0; i < flags.size(); ++i)
431  if (pointMarker == flags[i])
432  {
433  ++check;
434  }
435  return (check > 0) ? true : false;
436 }
437 
438 
439 template <typename mesh_Type>
440 double RBFrescaledScalar<mesh_Type>::rbf (double x1, double y1, double z1, double x2, double y2, double z2, double radius)
441 {
442  double distance = sqrt ( pow (x1 - x2, 2) + pow (y1 - y2, 2) + pow (z1 - z2, 2) );
443  return pow (1 - distance / radius, 4) * (4 * distance / radius + 1);
444 }
445 
446 template <typename mesh_Type>
448 {
449  *M_RhsF *= 0;
450  *M_RhsF = *newRhs;
451 }
452 
453 template <typename mesh_Type>
454 void RBFrescaledScalar<mesh_Type>::solution (vectorPtr_Type& Solution)
455 {
456  Solution = M_unknownField;
457 }
458 
459 
460 template <typename mesh_Type>
461 void RBFrescaledScalar<mesh_Type>::solutionrbf (vectorPtr_Type& Solution_rbf)
462 {
463  Solution_rbf = M_unknownField_rbf;
464 }
465 
466 //! Factory create function
467 template <typename mesh_Type>
469 {
470  return new RBFrescaledScalar< mesh_Type > ();
471 }
472 namespace
473 {
476 }
477 
478 } // Namespace LifeV
479 
480 #endif /* RBFRESCALEDSCALAR_H */
std::shared_ptr< MapEpetra > mapPtr_Type
GhostHandler< mesh_Type > neighbors_Type
void start()
Start the timer.
Definition: LifeChrono.hpp:93
void solution(vectorPtr_Type &Solution)
void setupRBFData(vectorPtr_Type KnownField, vectorPtr_Type UnknownField, GetPot datafile, parameterList_Type belosList)
idContainer_Type M_GIdsUnknownMesh
idContainer_Type M_GIdsKnownMesh
std::shared_ptr< vector_Type > vectorPtr_Type
std::shared_ptr< mesh_Type > meshPtr_Type
std::shared_ptr< prec_Type > precPtr_Type
LifeV::PreconditionerIfpack prec_Type
std::vector< int > flagContainer_Type
matrixPtr_Type M_interpolationOperator
Epetra_Import const & importer()
Getter for the Epetra_Import.
Definition: MapEpetra.cpp:394
void setRadius(double radius)
#define ASSERT(X, A)
Definition: LifeAssert.hpp:90
uint32_type ID
IDs.
Definition: LifeV.hpp:194
bool isInside(ID pointMarker, flagContainer_Type Flags)
std::unordered_set< ID > idContainer_Type
RBFInterpolation< mesh_Type > * createRBFrescaledScalar()
Factory create function.
std::shared_ptr< neighbors_Type > neighborsPtr_Type
Teuchos::RCP< Teuchos::ParameterList > parameterList_Type
std::shared_ptr< matrix_Type > matrixPtr_Type
double rbf(double x1, double y1, double z1, double x2, double y2, double z2, double radius)
void updateRhs(vectorPtr_Type newRhs)
std::shared_ptr< basePrec_Type > basePrecPtr_Type
void solutionrbf(vectorPtr_Type &Solution_rbf)
parameterList_Type M_belosList
void identifyNodes(meshPtr_Type LocalMesh, std::unordered_set< ID > &GID_nodes, vectorPtr_Type CheckVector)
void stop()
Stop the timer.
Definition: LifeChrono.hpp:100
void setup(meshPtr_Type fullMeshKnown, meshPtr_Type localMeshKnown, meshPtr_Type fullMeshUnknown, meshPtr_Type localMeshUnknown, flagContainer_Type flags)
LifeV::Preconditioner basePrec_Type
void importFromHDF5(std::string const &fileName="ghostmap")
Import neighbor lists to an hdf5 file.
MatrixEpetra< double > matrix_Type