Commit 9436807e authored by quinodoz's avatar quinodoz
Browse files

Changing the LocalDofPattern class:

- changing from int* to std::vector
- no more exposed members
- added accessors
- added constructors/destructors
- made 1D-2D-3D compliant
- corrected a bug for the P1isoP2
- added documentation

Changes in the other files are only to change from direct access to use of the accessors (same names, but with () )
parent 54165f20
......@@ -395,28 +395,28 @@ EpetraMap::setUp(const RefFE& refFE,
{
int indexBase = 1;
if (refFE.nbDofPerVertex)
if (refFE.nbDofPerVertex())
{
int numNode = repeatedNodeVector.size();
EpetraMap repeatedNodeMap(-1, numNode, &repeatedNodeVector[0], indexBase, _comm);
operator+=(repeatedNodeMap);
}
if (refFE.nbDofPerEdge)
if (refFE.nbDofPerEdge())
{
int numEdge = repeatedEdgeVector.size();
EpetraMap repeatedEdgeMap(-1, numEdge, &repeatedEdgeVector[0], indexBase, _comm);
operator+=(repeatedEdgeMap);
}
if (refFE.nbDofPerFace)
if (refFE.nbDofPerFace())
{
int numFace = repeatedFaceVector.size();
EpetraMap repeatedFaceMap(-1, numFace, &repeatedFaceVector[0], indexBase, _comm);
operator+=(repeatedFaceMap);
}
if (refFE.nbDofPerVolume)
if (refFE.nbDofPerVolume())
{
int numElem = repeatedVolumeVector.size();
EpetraMap repeatedElemMap(-1, numElem, &repeatedVolumeVector[0], indexBase, _comm);
......
......@@ -279,14 +279,14 @@ EpetraMap(const RefFE& refFE,
std::vector<int> repeatedFaceVector;
std::vector<int> repeatedVolumeVector;
if (refFE.nbDofPerVertex)
if (refFE.nbDofPerVertex())
{
repeatedNodeVector.reserve(mesh.numPoints());
for ( UInt ii = 1; ii <= mesh.numPoints(); ii++ )
repeatedNodeVector.push_back(mesh.pointList( ii ).id());
}
if (refFE.nbDofPerEdge)
if (refFE.nbDofPerEdge())
{
repeatedEdgeVector.reserve(mesh.numEdges());
......@@ -294,7 +294,7 @@ EpetraMap(const RefFE& refFE,
repeatedEdgeVector.push_back(mesh.edgeList( ii ).id());
}
if (refFE.nbDofPerFace)
if (refFE.nbDofPerFace())
{
repeatedFaceVector.reserve(mesh.numFaces());
......@@ -302,7 +302,7 @@ EpetraMap(const RefFE& refFE,
repeatedFaceVector.push_back(mesh.faceList( ii ).id());
}
if (refFE.nbDofPerVolume)
if (refFE.nbDofPerVolume())
{
repeatedVolumeVector.reserve(mesh.numVolumes());
......
......@@ -1552,7 +1552,7 @@ assemb_vec( EpetraVector& V, ElemVec& elvec, const LocalDofPattern& fe, const DO
// std::cout << "in assemb_vec" << std::endl;
UInt ig;
// UInt eleId = feId; //simplify.
for ( i = 0 ; i < fe.nbLocalDof ; i++ )
for ( i = 0 ; i < fe.nbLocalDof() ; i++ )
{ //! instead of CurrentFE::nbNode
ig = dof.localToGlobal( feId, i + 1 ) - 1 + iblock * totdof;
// std::cout << "i= " << i << std::endl;
......@@ -1583,7 +1583,7 @@ extract_vec( EpetraVector& V, ElemVec& elvec, const LocalDofPattern& fe, const D
int i;
// std::cout << "in assemb_vec" << std::endl;
UInt ig;
for ( i = 0 ; i < fe.nbLocalDof ; i++ )
for ( i = 0 ; i < fe.nbLocalDof() ; i++ )
{
ig = dof.localToGlobal( feId, i + 1 ) - 1 + iblock * totdof;
// std::cout << "i= " << i << std::endl;
......
......@@ -26,7 +26,7 @@ Dof::Dof( const LocalDofPattern& _fe, UInt off ) : fe( _fe ), _offset( off ), _t
_numFaces(0),_ltgByFace(),_gtlByFace()
{
//Getting the face
switch( _fe.nbLocalDof )
switch( _fe.nbLocalDof() )
{
case 4:
_fToP = LinearTetra::fToP;
......
......@@ -115,7 +115,7 @@ public:
//! The number of local Dof (nodes) in the finite element
inline UInt numLocalDof() const
{
return fe.nbLocalDof;
return fe.nbLocalDof();
}
//! Return the specified entries of the localToGlobal table
......@@ -226,7 +226,7 @@ Dof::Dof( Mesh& mesh, const LocalDofPattern& _fe, UInt off ) :
_gtlByFace()
{
//Getting the face
switch( _fe.nbLocalDof )
switch( _fe.nbLocalDof() )
{
case 4:
_fToP = LinearTetra::fToP;
......@@ -268,10 +268,10 @@ void Dof::update( Mesh& M )
typedef typename Mesh::ElementShape GeoShape;
// Some useful local variables, to save some typing
UInt nldpe = fe.nbDofPerEdge;
UInt nldpv = fe.nbDofPerVertex;
UInt nldpf = fe.nbDofPerFace;
UInt nldpV = fe.nbDofPerVolume;
UInt nldpe = fe.nbDofPerEdge();
UInt nldpv = fe.nbDofPerVertex();
UInt nldpf = fe.nbDofPerFace();
UInt nldpV = fe.nbDofPerVolume();
nlv = GeoShape::numVertices;
nle = GeoShape::numEdges;
......@@ -294,7 +294,7 @@ void Dof::update( Mesh& M )
UInt nldof = nldpV + nldpe * nle + nldpv * nlv + nldpf * nlf;
ASSERT_PRE( nldof == UInt( fe.nbLocalDof ), "Something wrong in FE specification" ) ;
ASSERT_PRE( nldof == UInt( fe.nbLocalDof() ), "Something wrong in FE specification" ) ;
_totalDof = nV * nldpV + ne * nldpe + nv * nldpv + nf * nldpf;
......
......@@ -70,7 +70,7 @@ public:
inline UInt numTotalDof() const {return _totalDof;}
//! The number of local Dof (nodes) in the finite element
inline UInt numLocalDof() const {return fe.nbLocalDof;}
inline UInt numLocalDof() const {return fe.nbLocalDof();}
//! Return the specified entries of the localToGlobal table
/*!
......
/*
This file is part of the LifeV library
Copyright (C) 2001,2002,2003,2004 EPFL, INRIA and Politecnico di Milano
......@@ -304,13 +305,13 @@ void DofInterface3Dto3D::_updateDofConnections( const Mesh& mesh1, const Dof& do
UInt nFaceV = GeoBShape::numVertices; // Number of face's vertices
UInt nFaceE = GeoBShape::numEdges; // Number of face's edges
UInt nDofpV1 = _refFE1->nbDofPerVertex; // number of Dof per vertices on mesh1
UInt nDofpE1 = _refFE1->nbDofPerEdge; // number of Dof per edges on mesh1
UInt nDofpF1 = _refFE1->nbDofPerFace; // number of Dof per faces on mesh1
UInt nDofpV1 = _refFE1->nbDofPerVertex(); // number of Dof per vertices on mesh1
UInt nDofpE1 = _refFE1->nbDofPerEdge(); // number of Dof per edges on mesh1
UInt nDofpF1 = _refFE1->nbDofPerFace(); // number of Dof per faces on mesh1
UInt nDofpV2 = _refFE2->nbDofPerVertex; // number of Dof per vertices on mesh2
UInt nDofpE2 = _refFE2->nbDofPerEdge; // number of Dof per edges on mesh2
UInt nDofpF2 = _refFE2->nbDofPerFace; // number of Dof per faces on mesh2
UInt nDofpV2 = _refFE2->nbDofPerVertex(); // number of Dof per vertices on mesh2
UInt nDofpE2 = _refFE2->nbDofPerEdge(); // number of Dof per edges on mesh2
UInt nDofpF2 = _refFE2->nbDofPerFace(); // number of Dof per faces on mesh2
UInt nElemV = GeoShape::numVertices; // Number of element's vertices
UInt nElemE = GeoShape::numEdges; // Number of element's edges
......
//@HEADER
/*
This file is part of the LifeV library
Copyright (C) 2001,2002,2003,2004 EPFL, INRIA and Politecnico di Milano
************************************************************************
This library 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 2.1 of the License, or (at your option) any later version.
This file is part of the LifeV Applications.
Copyright (C) 2001-2010 EPFL, Politecnico di Milano, INRIA
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
This library 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 2.1 of the
License, or (at your option) any later version.
This library 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 this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
USA
************************************************************************
*/
//@HEADER
/*!
@file
@brief Implementation of the localDofPattern class.
*/
#include <life/lifefem/localDofPattern.hpp>
namespace LifeV
{
LocalDofPattern::LocalDofPattern( int _nbLocalDof, int _nbDofPerVertex, int _nbDofPerEdge, int _nbDofPerFace,
int _nbDofPerVolume, PatternType _patternType ) :
nbLocalDof( _nbLocalDof ), nbDofPerVertex( _nbDofPerVertex ), nbDofPerEdge( _nbDofPerEdge ),
nbDofPerFace( _nbDofPerFace ), nbDofPerVolume( _nbDofPerVolume ),
patternType( _patternType )
// nbdof is total number of dof (not used for the initialization of a member)
// This is the 3D constructor
LocalDofPattern::LocalDofPattern( const UInt& _nbLocalDof, const UInt& _nbDofPerVertex,
const UInt& _nbDofPerEdge, const UInt& _nbDofPerFace,
const UInt& _nbDofPerVolume, const DofPatternType& _patternType ) :
M_dim(3), M_nbLocalDof( _nbLocalDof ), M_nbDofPerDimEntity(std::vector< UInt> (4)),
M_patternType( _patternType )
{
// Pattern
switch ( patternType )
// Store the location of the dofs
M_nbDofPerDimEntity[0]=_nbDofPerVertex;
M_nbDofPerDimEntity[1]=_nbDofPerEdge;
M_nbDofPerDimEntity[2]=_nbDofPerFace;
M_nbDofPerDimEntity[3]=_nbDofPerVolume;
// Decide the pattern depending on the type
switch ( M_patternType )
{
case STANDARD_PATTERN:
{
// full element matrix
_nbPattern = nbLocalDof * nbLocalDof;
_nbDiag = nbLocalDof;
_nbUpper = nbLocalDof * ( nbLocalDof - 1 ) / 2;
_patternFirst = new int[ _nbPattern ];
_patternSecond = new int[ _nbPattern ];
for ( int i = 0;i < nbLocalDof;i++ )
_patternFirst[ i ] = _patternSecond[ i ] = i;
int ip = nbLocalDof;
for ( int i = 0;i < nbLocalDof - 1;i++ )
{
for ( int j = i + 1;j < nbLocalDof;j++ )
{
_patternFirst[ ip ] = i;
_patternSecond[ ip ] = j;
ip++;
}
}
for ( int i = 1;i < nbLocalDof;i++ )
{
for ( int j = 0;j < i;j++ )
{
_patternFirst[ ip ] = i;
_patternSecond[ ip ] = j;
ip++;
}
}
break;
}
{
setupStandardPattern();
break;
}
default:
{
std::ostringstream _err_msg;
_err_msg << "Pattern " << M_patternType << " not available for " << M_dim << "D. ";
ERROR_MSG( _err_msg.str().c_str() );
};
};
}
// This is the 2D constructor
LocalDofPattern::LocalDofPattern( const UInt& _nbLocalDof, const UInt& _nbDofPerVertex,
const UInt& _nbDofPerEdge, const UInt& _nbDofPerFace,
const DofPatternType& _patternType ) :
M_dim(2), M_nbLocalDof( _nbLocalDof ), M_nbDofPerDimEntity(std::vector< UInt> (3)),
M_patternType( _patternType )
{
// Store the location of the dofs
M_nbDofPerDimEntity[0]=_nbDofPerVertex;
M_nbDofPerDimEntity[1]=_nbDofPerEdge;
M_nbDofPerDimEntity[2]=_nbDofPerFace;
// Decide the pattern depending on the type
switch ( M_patternType )
{
case STANDARD_PATTERN:
{
setupStandardPattern();
break;
}
case P1ISOP2_TRIA_PATTERN:
{
_nbPattern = 24;
_nbDiag = 6;
_nbUpper = 9;
typedef std::pair<int, int> COO;
COO pattern_p1isop2_tria[ 24 ] = {
COO( 0, 0 ), COO( 1, 1 ), COO( 2, 2 ), COO( 3, 3 ), COO( 4, 4 ), COO( 5, 5 ), // diagonal entries
COO( 0, 3 ), COO( 0, 5 ), COO( 1, 3 ), COO( 1, 4 ), COO( 2, 4 ), COO( 2, 5 ), COO( 3, 4 ), COO( 3, 5 ), COO( 4, 5 ), //upper entries
COO( 3, 0 ), COO( 3, 1 ), COO( 4, 1 ), COO( 4, 2 ), COO( 4, 3 ), COO( 5, 0 ), COO( 5, 2 ), COO( 5, 3 ), COO( 5, 4 ) // lower entries
};
for ( int i = 0;i < 24;i++ )
{
_patternFirst[ i ] = pattern_p1isop2_tria[ i ].first;
_patternSecond[ i ] = pattern_p1isop2_tria[ i ].second;
}
}
{
setupP1isoP2TriaPattern();
break;
}
default:
std::ostringstream _err_msg;
_err_msg << "Unknown pattern " << patternType << "I cannot build local pattern!";
{
std::ostringstream _err_msg;
_err_msg << "Pattern " << M_patternType << " not available for " << M_dim << "D. ";
ERROR_MSG( _err_msg.str().c_str() );
};
};
}
// This is the 1D constructor
LocalDofPattern::LocalDofPattern( const UInt& _nbLocalDof, const UInt& _nbDofPerVertex,
const UInt& _nbDofPerEdge, const DofPatternType& _patternType ) :
M_dim(1), M_nbLocalDof( _nbLocalDof ), M_nbDofPerDimEntity(std::vector< UInt> (2)),
M_patternType( _patternType )
{
// Store the location of the dofs
M_nbDofPerDimEntity[0]=_nbDofPerVertex;
M_nbDofPerDimEntity[1]=_nbDofPerEdge;
// Decide the pattern depending on the type
switch ( M_patternType )
{
case STANDARD_PATTERN:
{
setupStandardPattern();
break;
}
case P1ISOP2_SEG_PATTERN:
{
setupP1isoP2SegPattern();
break;
}
default:
{
std::ostringstream _err_msg;
_err_msg << "Pattern " << M_patternType << " not available for " << M_dim << "D. ";
ERROR_MSG( _err_msg.str().c_str() );
};
}; // end of the switch
};
// The copy constructor
LocalDofPattern::LocalDofPattern( const LocalDofPattern& _localDofPattern) :
M_dim(_localDofPattern.M_dim),
M_nbLocalDof (_localDofPattern.M_nbLocalDof ),
M_nbDofPerDimEntity (_localDofPattern.M_nbDofPerDimEntity),
M_patternType (_localDofPattern.M_patternType ),
M_pattern (_localDofPattern.M_pattern),
M_nbPattern(_localDofPattern.M_nbPattern),
M_nbDiag (_localDofPattern.M_nbDiag ),
M_nbUpper (_localDofPattern.M_nbUpper )
{};
void LocalDofPattern::showMe( std::ostream& output) const
{
output << " Size of the pattern : " << M_nbPattern << std::endl;
output << " Diag: " << M_nbDiag << " Upper: " << M_nbUpper << std::endl;
output << " Pattern type: " << M_patternType << std::endl;
for (UInt iter(0); iter< M_nbPattern; ++iter)
{
output << iter << " : " << M_pattern[iter].first << " - " << M_pattern[iter].second << std::endl;
};
}
void LocalDofPattern::setupStandardPattern()
{
// This is the standard pattern with all the
// degrees of freedom coupled together.
// Number of couplings
M_nbPattern = M_nbLocalDof * M_nbLocalDof;
M_nbDiag = M_nbLocalDof;
M_nbUpper = M_nbLocalDof * ( M_nbLocalDof - 1 ) / 2;
// initialization of the pattern
M_pattern = std::vector< std::pair< UInt, UInt > > (M_nbPattern);
// First, put the diagonal entries
for ( UInt i = 0;i < M_nbLocalDof;i++ )
{
M_pattern[i] = std::pair<UInt,UInt> (i,i);
}
// Upper diagonal entries
int ip ( M_nbLocalDof );
for ( UInt i = 0;i < M_nbLocalDof - 1;i++ )
{
for ( UInt j = i + 1;j < M_nbLocalDof;j++ )
{
M_pattern[ip] = std::pair<UInt,UInt> (i,j);
ip++;
}
}
// Lower diagonal entries
for ( UInt i = 1;i < M_nbLocalDof;i++ )
{
for ( UInt j = 0;j < i;j++ )
{
M_pattern[ip] = std::pair<UInt,UInt> (i,j);
ip++;
}
};
};
void LocalDofPattern::setupP1isoP2SegPattern()
{
// Some check to ensure consistency
ASSERT(M_nbDofPerDimEntity[0] == 1, " Inconsistent P1 iso P2 (Vertices)");
ASSERT(M_nbDofPerDimEntity[1] == 1, " Inconsistent P1 iso P2 (Edges)");
ASSERT(M_nbDofPerDimEntity[2] == 0, " Inconsistent P1 iso P2 (Faces)");
M_nbPattern = 7;
M_nbDiag = 3;
M_nbUpper = 2;
M_pattern = std::vector< std::pair < UInt,UInt > > (M_nbPattern);
// Diagonal entries
for (UInt diag(0); diag<3; ++diag)
{
M_pattern[diag] = std::pair<UInt,UInt> (diag,diag);
};
// Upper diagonal entries
M_pattern[3] = std::pair<UInt,UInt> (0,2);
M_pattern[4] = std::pair<UInt,UInt> (1,2);
// Lower diagonal entries
M_pattern[5] = std::pair<UInt,UInt> (2,0);
M_pattern[6] = std::pair<UInt,UInt> (2,1);
}
void LocalDofPattern::setupP1isoP2TriaPattern()
{
// Some check to ensure consistency
ASSERT(M_nbDofPerDimEntity[0] == 1, " Inconsistent P1 iso P2 (Vertices)");
ASSERT(M_nbDofPerDimEntity[1] == 1, " Inconsistent P1 iso P2 (Edges)");
ASSERT(M_nbDofPerDimEntity[2] == 0, " Inconsistent P1 iso P2 (Faces)");
M_nbPattern = 24;
M_nbDiag = 6;
M_nbUpper = 9;
M_pattern = std::vector< std::pair < UInt,UInt > > (M_nbPattern);
// Diagonal entries
for (UInt diag(0); diag<6; ++diag)
{
M_pattern[diag] = std::pair<UInt,UInt> (diag,diag);
};
// Upper diagonal entries
M_pattern[6] = std::pair<UInt,UInt> (0,3);
M_pattern[7] = std::pair<UInt,UInt> (0,5);
M_pattern[8] = std::pair<UInt,UInt> (1,3);
M_pattern[9] = std::pair<UInt,UInt> (1,4);
M_pattern[10] = std::pair<UInt,UInt> (2,4);
M_pattern[11] = std::pair<UInt,UInt> (2,5);
M_pattern[12] = std::pair<UInt,UInt> (3,4);
M_pattern[13] = std::pair<UInt,UInt> (3,5);
M_pattern[14] = std::pair<UInt,UInt> (4,5);
// Lower diagonal entries
M_pattern[15] = std::pair<UInt,UInt> (3,0);
M_pattern[16] = std::pair<UInt,UInt> (3,1);
M_pattern[17] = std::pair<UInt,UInt> (4,1);
M_pattern[18] = std::pair<UInt,UInt> (4,2);
M_pattern[19] = std::pair<UInt,UInt> (4,3);
M_pattern[20] = std::pair<UInt,UInt> (5,0);
M_pattern[21] = std::pair<UInt,UInt> (5,2);
M_pattern[22] = std::pair<UInt,UInt> (5,3);
M_pattern[23] = std::pair<UInt,UInt> (5,4);
}
} // end of the namespace LifeV
//@HEADER
/*
This file is part of the LifeV library
Copyright (C) 2001,2002,2003,2004 EPFL, INRIA and Politecnico di Milano
************************************************************************
This library 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 2.1 of the License, or (at your option) any later version.
This file is part of the LifeV Applications.
Copyright (C) 2001-2010 EPFL, Politecnico di Milano, INRIA
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
This library 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 2.1 of the
License, or (at your option) any later version.
This library 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 this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
USA
************************************************************************
*/
//@HEADER
/*!
@file
@brief This file contains the definition of the LocalDofPattern class.
*/
#ifndef _LOCAL_DOF_PATTERN_HH
#define _LOCAL_DOF_PATTERN_HH
......@@ -24,56 +38,252 @@
namespace LifeV
{
/*! Local pattern type
//! Local pattern type
/*!
This enum allows to distinguish the normal standard local pattern, which is a full pattern involving all
degrees of freedom to special patterns. It is stared in LocalDofPattern for later use by Dof
degrees of freedom to special patterns. It is stored in LocalDofPattern for later use by Dof
*/
enum PatternType {STANDARD_PATTERN = 1, P1ISOP2_TRIA_PATTERN = 2};
enum DofPatternType {STANDARD_PATTERN = 1,
P1ISOP2_SEG_PATTERN = 2,
P1ISOP2_TRIA_PATTERN = 3};
//! LocalDofPattern - A class to store the "couplings" between the basis functions
/*!
The aim of this class is to store the way the basis functions couple one with each other. This might seem useless,
however, some "advanced" finite elements require this structure.
For example, consider the P1-iso-P2 element in 2D. This finite element is composed of 6 basis functions, based on the
nodes with the same numerotation as the P2 element. The reference triangle is split into 4 subtriangles using the
nodes on the faces. Each basis function is build such that it is 1 on its node, 0 on all the other nodes and such that
it is linear in each subtriangle.
@see in "Numerical Approximation of Partial Differential Equations" by A. Quarteroni and A. Valli, p.311, for an
illustration and further informations.
With this definition of the P1-iso-P2 finite element, we see that the basis functions 1 and 2 have no common support.
So, they are not directly coupled.
In order to represent the couplings between the basis functions, we can use a matrix \f$ C \f$ such that
\f$ C_{ij} = 1 \f$ if the basis functions \f$ i \f$ and \f$ j \f$ have a common support, otherwise it is \f$ 0 \f$.
This matrix is symetric.
For the P1-iso-P2 element, this matrix would be:
\f[