LifeV
test_interpolate.cpp
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 Interpolate test
30 
31  @author Mauro Perego <mperego@fsu.edu>
32  @contributor
33  @maintainer Mauro Perego <mperego@fsu.edu>
34 
35  @date 06-30-2010
36 
37 The program tests the interpolation methods between different finite elements (mostly between scalar continuous finite elements).
38 Also it test the interpolation of an analytical function into a finite element space.
39  */
40 
41 
42 #include <Epetra_ConfigDefs.h>
43 #ifdef EPETRA_MPI
44 #include <mpi.h>
45 #include <Epetra_MpiComm.h>
46 #else
47 #include <Epetra_SerialComm.h>
48 #endif
49 
50 
51 
52 #include <lifev/core/fem/ReferenceFE.hpp>
53 #include <lifev/core/fem/QuadratureRule.hpp>
54 #include <lifev/core/fem/CurrentFE.hpp>
55 #include <lifev/core/filter/GetPot.hpp>
56 #include <lifev/core/mesh/RegionMesh.hpp>
57 #include <lifev/core/mesh/RegionMesh3DStructured.hpp>
58 #include <lifev/core/fem/DOF.hpp>
59 #include <lifev/core/filter/MeshWriter.hpp>
60 
61 #include <lifev/core/mesh/MeshData.hpp>
62 
63 
64 #include "test_interpolate.hpp"
65 
66 using namespace LifeV;
67 
68 int main (int argc, char** argv )
69 {
70  typedef FESpace < RegionMesh<LinearTetra>, MapEpetra > FESpaceTetra_Type;
71  typedef std::shared_ptr < FESpaceTetra_Type > FESpaceTetraPtr_Type;
72 
73  typedef FESpace < RegionMesh<LinearHexa>, MapEpetra > FESpaceHexa_Type;
74  typedef std::shared_ptr < FESpaceHexa_Type > FESpaceHexaPtr_Type;
75 
76 
77 #ifdef HAVE_MPI
78  MPI_Init (&argc, &argv);
79  std::shared_ptr<Epetra_Comm> Comm (new Epetra_MpiComm (MPI_COMM_WORLD) );
80 #else
81  std::shared_ptr<Epetra_Comm> Comm (new Epetra_SerialComm);
82 #endif
83 
84  UInt verbose = (Comm->MyPID() == 0);
85 
86  bool check (true);
87 
88 
89  //definition of Array of Errors which will be used to check the correctness of the interpolate test.
90 
91  const Real errArrayBilinear[2] = { 0.0312819802, 0. };
92 
93  const Real errArrayQuadratic[12] = { 0.0136247667, 0.0005088372, 0.0005577494, 0.0005088372,
94  0.0136172446, 0.0005088372, 0.0004270717, 0.0005088372,
95  0.0136172446, 0.0005088372, 0.0004270717, 0.
96  };
97 
98  const Real errArrayBubble[12] = { 0.0094702745, 3.584186e-10, 3.67611e-10, 3.584186e-10,
99  0.0094702745, 3.584186e-10, 0., 3.584186e-10,
100  0.0094702745, 3.584186e-10, 3.67611e-10, 3.584186e-10
101  };
102 
103  const Real errArrayLinear[12] = { 0.010437463587, 0., 0., 0.,
104  0.010437463587, 0., 0., 0.,
105  0.010437463587, 0., 0., 0.
106  };
107 
108 
109  const std::string stringArrayP[12] = { "P1 -> P0 ", "P1 -> P1 ", "P1 -> P1b", "P1 -> P2 ",
110  "P1b -> P0 ", "P1b -> P1 ", "P1b -> P1b", "P1b -> P2 ",
111  "P2 -> P0 ", "P2 -> P1 ", "P2 -> P1b", "P2 -> P2 "
112  };
113 
114  const std::string stringArrayQ[2] = {"Q1 -> Q0 ", "Q1 -> Q1 "};
115 
116 
117  // Import/Generate an hexahedral and a Tetrahedral mesh.
118 
119  std::shared_ptr<RegionMesh<LinearTetra> > fullMeshTetraPtr ( new RegionMesh<LinearTetra> ( Comm ) );
120  std::shared_ptr<RegionMesh<LinearHexa> > fullMeshHexaPtr ( new RegionMesh<LinearHexa> ( Comm ) );
121  UInt nEl (10);
122  GetPot dataFile ("./data");
123  MeshData meshData (dataFile, "interpolate/space_discretization");
124  readMesh (*fullMeshHexaPtr, meshData);
125  regularMesh3D ( *fullMeshTetraPtr, 1, nEl, nEl, nEl);
126 
127 
128  // Partition the meshes using ParMetis
129  std::shared_ptr<RegionMesh<LinearHexa> > localMeshHexaPtr;
130  {
131  // Create the partitioner
132  MeshPartitioner< RegionMesh<LinearHexa> > meshPartHexa;
133 
134  // Partition the mesh using ParMetis
135  meshPartHexa.doPartition ( fullMeshHexaPtr, Comm );
136 
137  // Get the local mesh
138  localMeshHexaPtr = meshPartHexa.meshPartition();
139  }
140  std::shared_ptr<RegionMesh<LinearTetra> > localMeshTetraPtr;
141  {
142  // Create the partitioner
143  MeshPartitioner< RegionMesh<LinearTetra> > meshPartTetra;
144 
145  // Partition the mesh using ParMetis
146  meshPartTetra.doPartition ( fullMeshTetraPtr, Comm );
147 
148  // Get the local mesh
149  localMeshTetraPtr = meshPartTetra.meshPartition();
150  }
151 
152  //Building finite element spaces
153 
154  //Finite element space of the first scalar field - P0
155  FESpaceTetraPtr_Type feSpaceP0 ( new FESpaceTetra_Type ( localMeshTetraPtr, feTetraP0, quadRuleTetra4pt,
156  quadRuleTria4pt, 3, Comm ) );
157 
158  //Finite element space of the second scalar field - P1
159  FESpaceTetraPtr_Type feSpaceP1 ( new FESpaceTetra_Type ( localMeshTetraPtr, feTetraP1, quadRuleTetra4pt,
160  quadRuleTria4pt, 3, Comm ) );
161 
162  // Finite element space of the second scalar field - P1bubble
163  FESpaceTetraPtr_Type feSpaceP1Bubble ( new FESpaceTetra_Type ( localMeshTetraPtr, feTetraP1bubble, quadRuleTetra15pt,
164  quadRuleTria4pt, 3, Comm ) );
165 
166  // Finite element space of the second scalar field - P2
167  FESpaceTetraPtr_Type feSpaceP2 ( new FESpaceTetra_Type ( localMeshTetraPtr, feTetraP2, quadRuleTetra4pt,
168  quadRuleTria4pt, 3, Comm ) );
169 
170  // Finite element space of the first scalar field - Q0
171  FESpaceHexaPtr_Type feSpaceQ0 ( new FESpaceHexa_Type ( localMeshHexaPtr, feHexaQ0, quadRuleHexa8pt,
172  quadRuleQuad4pt, 2, Comm ) );
173 
174  // Finite element space of the second scalar field - Q1
175  FESpaceHexaPtr_Type feSpaceQ1 ( new FESpaceHexa_Type ( localMeshHexaPtr, feHexaQ1, quadRuleHexa8pt,
176  quadRuleQuad4pt, 2, Comm ) );
177 
178  //vectors containing the original and final FE spaces
179  //(FE vectors will be interpolated from original FE spaces into final FE Spaces)
180  std::vector<FESpaceHexaPtr_Type > originalFeSpaceHexaVec (1), finalFeSpaceHexaVec (2);
181  std::vector<FESpaceTetraPtr_Type> originalFeSpaceTetraVec (3), finalFeSpaceTetraVec (4);
182 
183  originalFeSpaceHexaVec[0] = feSpaceQ1;
184  finalFeSpaceHexaVec[0] = feSpaceQ0;
185  finalFeSpaceHexaVec[1] = feSpaceQ1;
186 
187  originalFeSpaceTetraVec[0] = feSpaceP1;
188  originalFeSpaceTetraVec[1] = feSpaceP1Bubble;
189  originalFeSpaceTetraVec[2] = feSpaceP2;
190 
191  finalFeSpaceTetraVec[0] = feSpaceP0;
192  finalFeSpaceTetraVec[1] = feSpaceP1;
193  finalFeSpaceTetraVec[2] = feSpaceP1Bubble;
194  finalFeSpaceTetraVec[3] = feSpaceP2;
195 
196  Real time = 0.5;
197 
198  if (verbose)
199  {
200  std::cout << "\nA bilinear function is interpolated into Q1 vector. \nThen this FE vector is interpolated into Q0 and Q1 vectors. \nThese are the errors with respect to the analytical solution.\n";
201  }
202  check = check_interpolate (originalFeSpaceHexaVec, finalFeSpaceHexaVec, Unique, bilinearFunction, errArrayBilinear, stringArrayQ, 1e-10, time, verbose);
203 
204 
205  if (verbose)
206  std::cout << "\nA linear function is interpolated into P1, P1b, P2 vectors. "
207  "\nThen these FE vectors are interpolated into the following finite elements \nand error with respect to the analytic solution are reported.\n";
208 
209  check &= check_interpolate (originalFeSpaceTetraVec, finalFeSpaceTetraVec, Repeated, linearFunction, errArrayLinear, stringArrayP, 1e-10, time, verbose);
210 
211 
212  if (verbose)
213  std::cout << "\nA quadratic function is interpolated into P1, P1b, P2 vectors. "
214  "\nThen these FE vectors are interpolated into the following finite elements \nand error with respect to the analytic solution are reported.\n";
215 
216  check &= check_interpolate (originalFeSpaceTetraVec, finalFeSpaceTetraVec, Unique, quadraticFunction, errArrayQuadratic, stringArrayP, 1e-10, time, verbose);
217 
218 
219  if (verbose)
220  std::cout << "\nA linear bubble function is interpolated into P1, P1b, P2 vectors. "
221  "\nThen these FE vectors are interpolated into the following finite elements \nand error with respect to the analytic solution are reported.\n";
222 
223  check &= check_interpolate (originalFeSpaceTetraVec, finalFeSpaceTetraVec, Repeated, linearBubbleFunction, errArrayBubble, stringArrayP, 1e-10, time, verbose);
224 
225 
226 #ifdef HAVE_MPI
227  std::cout << "MPI Finalization" << std::endl;
228  MPI_Finalize();
229 #endif
230 
231  if (verbose)
232  {
233  if (check)
234  {
235  std::cout << "\nTEST INTERPOLATE WAS SUCCESSFUL.\n\n";
236  }
237  else
238  {
239  std::cout << "\nTEST INTERPOLATE FAILED.\n\n";
240  }
241  }
242 
243  if (check)
244  {
245  return EXIT_SUCCESS;
246  }
247  else
248  {
249  return EXIT_FAILURE;
250  }
251 }//end main
252 
253 
254 Real linearFunction (const Real& t, const Real& x, const Real& y, const Real& z, const ID& ic)
255 {
256  return 2 * x + -z + 4 + y * (ic + 1) + t;
257 }
258 
259 //linear function containing a bubble. The bubble is defined on a particular mesh tetrahedra.
260 Real linearBubbleFunction (const Real& t, const Real& x, const Real& y, const Real& z, const ID& ic)
261 {
262  Real xx = std::max (x, 0.);
263  Real yy = std::max (y - x, 0.);
264  Real zz = std::max (z - y, 0.);
265  return 2 * x + -z + 4 + y * ic + t + xx * yy * zz * std::max (0.1 - xx - yy - zz, 0.);
266 }
267 
268 Real quadraticFunction (const Real& t, const Real& x, const Real& y, const Real& z, const ID& ic)
269 {
270  return 2 * x * x + t * y * y - z * z + x * y - x * z + 2 * y * z - x + y * (ic + 1) - z - 8;
271 }
272 
273 Real bilinearFunction (const Real& /*t*/, const Real& x, const Real& y, const Real& z, const ID& ic)
274 {
275  return x * y - x * z + 2 * y * z - x + y * (ic + 1) - z - 8;
276 }
main
int main(int argc, char **argv)
Definition: lifev/structure/testsuite/evaluateNodalETA/main.cpp:921