Commit ec9f30b1 authored by Luca Pasquale's avatar Luca Pasquale
Browse files

Created tutorial 9_ETA_gradient_interpolation

Tutorial 9_ETA_gradient_interpolation does the same things as the test
ETA_interpolate_gradient_2D, but in 3D
Fixed various comments in other ETA tests
parent 5edadbee
......@@ -182,8 +182,8 @@ ETA_ADR1DTest::run()
// ---------------------------------------------------------------
// We interpolate then the advection function of the mesh at hand.
// This is performed with the classical FESpace only.
// We then interpolate the advection function.
// This is can only be performed with the classical FESpace.
// ---------------------------------------------------------------
if (verbose)
......@@ -224,6 +224,7 @@ ETA_ADR1DTest::run()
// ---------------------------------------------------------------
// Definition of the RHS
// ---------------------------------------------------------------
if (verbose)
{
......
......@@ -186,8 +186,8 @@ ETA_ADR2DTest::run()
// ---------------------------------------------------------------
// We interpolate then the advection function of the mesh at hand.
// This is performed with the classical FESpace only.
// We then interpolate the advection function.
// This is can only be performed with the classical FESpace.
// ---------------------------------------------------------------
if (verbose)
......@@ -228,6 +228,7 @@ ETA_ADR2DTest::run()
// ---------------------------------------------------------------
// Definition of the RHS
// ---------------------------------------------------------------
if (verbose)
{
......
......@@ -59,7 +59,7 @@ typedef VectorEpetra vector_Type;
// ===================================================
// ---------------------------------------------------------------
// We define then a function whose gradient is
// We define a function whose gradient is
// (1 0)
// (0 2)
// ---------------------------------------------------------------
......@@ -168,13 +168,13 @@ ETA_InterpolateGradient2DTest::run()
// ---------------------------------------------------------------
// We interpolate then the advection function of the mesh at hand.
// This is performed with the classical FESpace only.
// We interpolate then the function.
// This can only be performed with the classical FESpace.
// ---------------------------------------------------------------
if (verbose)
{
std::cout << " -- Interpolating the solution field ... " << std::flush;
std::cout << " -- Interpolating the function ... " << std::flush;
}
vector_Type uInterpolated (uSpace->map(), Unique);
......@@ -188,8 +188,8 @@ ETA_InterpolateGradient2DTest::run()
// ---------------------------------------------------------------
// We build define the SmallVector used to extract the trace
// and variable used to store the result of the integration
// We build the MatrixSmall used to extract the trace
// and a variable used to store the result of the integration
// ---------------------------------------------------------------
MatrixSmall<2, 2> identityMatrix;
......@@ -243,10 +243,8 @@ ETA_InterpolateGradient2DTest::run()
std::cout << " Time : " << ETChrono.diff() << std::endl;
}
// ---------------------------------------------------------------
// We finally need to check that both yield the same matrix. In
// that aim, we need to finalize both matrices.
// Integrals computed on each processor must be summed together
// ---------------------------------------------------------------
if (verbose)
......@@ -265,8 +263,8 @@ ETA_InterpolateGradient2DTest::run()
}
// ---------------------------------------------------------------
// We now compute the error as the difference the integral
// computed and the exact value expected (3*4=12)
// We now compute the error as the difference between the integral
// computed and the exact value expected ( (1+2)*4 = 12 )
// ---------------------------------------------------------------
if (verbose)
......
INCLUDE(TribitsAddExecutableAndTest)
INCLUDE(TribitsCopyFilesToBinaryDir)
INCLUDE_DIRECTORIES(${CMAKE_SOURCE_DIR})
TRIBITS_ADD_EXECUTABLE_AND_TEST(
9_ETA_gradient_interpolation
SOURCES main.cpp
ARGS -c
NUM_MPI_PROCS 2
COMM serial mpi
)
//@HEADER
/*
*******************************************************************************
Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
This file is part of LifeV.
LifeV is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
LifeV is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with LifeV. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************
*/
//@HEADER
/*!
@file
@brief Tutorial for gradient interpolation.
@author Luca Pasquale <luca.pasquale@mail.polimi.it>
@date 2013-05
In this tutorial we use two expressions available when
integrating with Expression Templates:
- interpoaltion of the gradient of a given function
- multiplication by a MatrixSmall
We do this by multiplying the gradient by the identity matrix
(i.e. we're integrating the divergence)
Tutorials that should be read before: 1,3
*/
// ---------------------------------------------------------------
// We reuse the same files as in the previous tutorial. We add the
// chrono to measure the timings.
// ---------------------------------------------------------------
#pragma GCC diagnostic ignored "-Wunused-variable"
#pragma GCC diagnostic ignored "-Wunused-parameter"
#include <Epetra_ConfigDefs.h>
#ifdef EPETRA_MPI
#include <mpi.h>
#include <Epetra_MpiComm.h>
#else
#include <Epetra_SerialComm.h>
#endif
#pragma GCC diagnostic warning "-Wunused-variable"
#pragma GCC diagnostic warning "-Wunused-parameter"
#include <lifev/core/LifeV.hpp>
#include <lifev/core/mesh/MeshPartitioner.hpp>
#include <lifev/core/mesh/RegionMesh3DStructured.hpp>
#include <lifev/core/mesh/RegionMesh.hpp>
#include <lifev/core/array/MatrixEpetra.hpp>
#include <lifev/core/array/VectorEpetra.hpp>
#include <lifev/eta/fem/ETFESpace.hpp>
#include <lifev/eta/expression/Integrate.hpp>
#include <lifev/core/fem/FESpace.hpp>
#include <lifev/core/util/LifeChrono.hpp>
#include <boost/shared_ptr.hpp>
// ---------------------------------------------------------------
// We work in the LifeV namespace and define the mesh, matrix and
// vector types that we will need several times.
// ---------------------------------------------------------------
using namespace LifeV;
typedef RegionMesh<LinearTetra> mesh_Type;
typedef MatrixEpetra<Real> matrix_Type;
typedef VectorEpetra vector_Type;
// ===================================================
//! Functions
// ===================================================
// ---------------------------------------------------------------
// We define a function whose gradient is
// (1 0 0)
// (0 5 0)
// (0 0 2)
// ---------------------------------------------------------------
Real uFct ( const Real& /* t */, const Real& x , const Real& y , const Real& z , const ID& i )
{
if (i == 0)
{
return x;
}
else if (i == 1)
{
return 5 * y;
}
return 2 * z;
}
// ---------------------------------------------------------------
// As usual, we start by the MPI communicator, the definition of
// the mesh and its partitioning.
// ---------------------------------------------------------------
int main ( int argc, char** argv )
{
#ifdef HAVE_MPI
MPI_Init (&argc, &argv);
boost::shared_ptr<Epetra_Comm> Comm (new Epetra_MpiComm (MPI_COMM_WORLD) );
#else
boost::shared_ptr<Epetra_Comm> Comm (new Epetra_SerialComm);
#endif
const bool verbose (Comm->MyPID() == 0);
if (verbose)
{
std::cout << " -- Building and partitioning the mesh ... " << std::flush;
}
const UInt Nelements (10);
boost::shared_ptr< mesh_Type > fullMeshPtr (new mesh_Type);
regularMesh3D ( *fullMeshPtr, 1, Nelements, Nelements, Nelements, false,
2.0, 2.0, 2.0,
-1.0, -1.0, -1.0);
MeshPartitioner< mesh_Type > meshPart (fullMeshPtr, Comm);
boost::shared_ptr< mesh_Type > meshPtr (meshPart.meshPartition() );
fullMeshPtr.reset();
if (verbose)
{
std::cout << " done ! " << std::endl;
}
// ---------------------------------------------------------------
// We define then the ETFESpace, that we suppose scalar in this
// case. We also need a standard FESpace to perform the
// interpolation.
// ---------------------------------------------------------------
if (verbose)
{
std::cout << " -- Building FESpaces ... " << std::flush;
}
std::string uOrder ("P1");
boost::shared_ptr<FESpace< mesh_Type, MapEpetra > > uSpace
( new FESpace< mesh_Type, MapEpetra > (meshPtr, uOrder, 3, Comm) );
if (verbose)
{
std::cout << " done ! " << std::endl;
}
if (verbose)
{
std::cout << " ---> Dofs: " << uSpace->dof().numTotalDof() << std::endl;
}
if (verbose)
{
std::cout << " -- Building ETFESpaces ... " << std::flush;
}
boost::shared_ptr<ETFESpace< mesh_Type, MapEpetra, 3, 3 > > ETuSpace
( new ETFESpace< mesh_Type, MapEpetra, 3, 3 > (meshPtr, & (uSpace->refFE() ), & (uSpace->fe().geoMap() ), Comm) );
if (verbose)
{
std::cout << " done ! " << std::endl;
}
if (verbose)
{
std::cout << " ---> Dofs: " << ETuSpace->dof().numTotalDof() << std::endl;
}
// ---------------------------------------------------------------
// We interpolate then the function.
// This can only be performed with the classical FESpace.
// ---------------------------------------------------------------
if (verbose)
{
std::cout << " -- Interpolating the function ... " << std::flush;
}
vector_Type uInterpolated (uSpace->map(), Unique);
uSpace->interpolate (static_cast<FESpace< mesh_Type, MapEpetra >::function_Type> (uFct), uInterpolated, 0.0);
vector_Type uInterpolatedRepeated (uInterpolated, Repeated);
if (verbose)
{
std::cout << " done! " << std::endl;
}
// ---------------------------------------------------------------
// We build the MatrixSmall used to extract the trace
// and a variable used to store the result of the integration.
// We only need to specify the diagonal elements of the MatrixSmall
// because its default constructor sets all the components to 0
// ---------------------------------------------------------------
MatrixSmall<3, 3> identityMatrix;
identityMatrix (0, 0) = 1;
identityMatrix (1, 1) = 1;
identityMatrix (2, 2) = 1;
Real ETintegral (0);
if (verbose)
{
std::cout << " done! " << std::endl;
}
// ---------------------------------------------------------------
// We integrate on the domain the trace of the gradient
// (1+5+2 in this case)
// ---------------------------------------------------------------
LifeChrono ETChrono;
ETChrono.start();
if (verbose)
{
std::cout << " -- ET assembly ... " << std::flush;
}
{
using namespace ExpressionAssembly;
integrate ( elements (ETuSpace->mesh() ),
uSpace->qr(),
dot ( grad (ETuSpace, uInterpolatedRepeated) , value (identityMatrix) )
)
>> ETintegral;
}
ETChrono.stop();
if (verbose)
{
std::cout << " done! " << std::endl;
}
if (verbose)
{
std::cout << " Time : " << ETChrono.diff() << std::endl;
}
// ---------------------------------------------------------------
// Integrals computed on each processor must be summed together
// ---------------------------------------------------------------
if (verbose)
{
std::cout << " -- Broadcasting and summing integrals across processes ... " << std::flush;
}
Real globalIntegral (0.0);
Comm->Barrier();
Comm->SumAll (&ETintegral, &globalIntegral, 1);
if (verbose)
{
std::cout << " done ! " << std::endl;
}
// ---------------------------------------------------------------
// We now compute the error as the difference between the integral
// computed and the exact value expected ( (1+5+2)*8 = 64 )
// ---------------------------------------------------------------
if (verbose)
{
std::cout << " -- Computing the error ... " << std::flush;
}
Real error = std::abs (globalIntegral - 64.0);
if (verbose)
{
std::cout << "Error: " << error << std::endl;
}
if (verbose)
{
std::cout << " done ! " << std::endl;
}
#ifdef HAVE_MPI
MPI_Finalize();
#endif
Real tolerance (1e-10);
if (error < tolerance)
{
return ( EXIT_SUCCESS );
}
return ( EXIT_FAILURE );
}
......@@ -11,4 +11,5 @@ ADD_SUBDIRECTORIES(
6_ETA_functor
7_ETA_blocks
8_ETA_block_manip
9_ETA_gradient_interpolation
)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment