Commit 8941429c authored by passerini's avatar passerini
Browse files

lifev-mox is now merged into lifev-head

parent 741e5d58
......@@ -63,3 +63,22 @@ PRESENT_MERGE_mox_to_head_2007_11_16
LATEST_MERGE_head_to_mox_2007_11_16
LATEST_MERGE_mox_to_head_2007_11_16
PRESENT_MERGE_head_to_parallel_2008_06_12
PRESENT_MERGE_parallel_to_head_2008_06_12
LATEST_MERGE_head_to_parallel_2008_06_12
LATEST_MERGE_parallel_to_head_2008_06_12
PRESENT_MERGE_head_to_mox_2008_06_24
PRESENT_MERGE_mox_to_head_2008_06_24
LATEST_MERGE_head_to_mox_2008_06_24
LATEST_MERGE_mox_to_head_2008_06_24
PRESENT_MERGE_head_to_parallel_2008_07_02
PRESENT_MERGE_parallel_to_head_2008_07_02
LATEST_MERGE_parallel_to_head_2008_07_02
LATEST_MERGE_head_to_parallel_2008_07_02
PRESENT_MERGE_head_to_mox_2008_12_02
PRESENT_MERGE_mox_to_head_2008_12_02
LATEST_MERGE_head_to_mox_2008_12_02
LATEST_MERGE_mox_to_head_2008_12_02
......@@ -2,19 +2,19 @@
# DO NOT EDIT : generated automatically by update-headers
#
LIFECORE_HDRS=
LIFECORE_HDRS=life/lifecore/dataString.hpp life/lifecore/application.hpp life/lifecore/GetPot.hpp life/lifecore/lifeassert.hpp life/lifecore/policy.hpp life/lifecore/fortran_wrap.hpp life/lifecore/debug.hpp life/lifecore/chrono.hpp life/lifecore/singleton.hpp life/lifecore/switch.hpp life/lifecore/typeInfo.hpp life/lifecore/lifeversion.hpp life/lifecore/lifemacros.hpp life/lifecore/factory.hpp life/lifecore/life.hpp life/lifecore/about.hpp life/lifecore/SmartAssert.hpp life/lifecore/util_string.hpp
LIFEARRAY_HDRS=
LIFEARRAY_HDRS=life/lifearray/CSRMatrix.hpp life/lifearray/EpetraVector.hpp life/lifearray/boostmatrix.hpp life/lifearray/RNM_op.hpp life/lifearray/tab.hpp life/lifearray/connectivity.hpp life/lifearray/elemVec.hpp life/lifearray/VBRMatrix.hpp life/lifearray/EpetraMatrix.hpp life/lifearray/SimpleVect.hpp life/lifearray/elemMat.hpp life/lifearray/pattern.hpp life/lifearray/RNM_tpl.hpp life/lifearray/MixedMatrix.hpp life/lifearray/RNM_opc.hpp life/lifearray/blockMatrix.hpp life/lifearray/variables.hpp life/lifearray/MSRMatrix.hpp life/lifearray/vecUnknown.hpp life/lifearray/tridiagMatrix.hpp life/lifearray/RNM.hpp life/lifearray/sparseArray.hpp
LIFEALG_HDRS=
LIFEALG_HDRS=life/lifealg/picard.hpp life/lifealg/givens.hpp life/lifealg/newton.hpp life/lifealg/SolverAztec.hpp life/lifealg/saddlePointCG.hpp life/lifealg/iteration.hpp life/lifealg/gmres.hpp life/lifealg/EpetraPreconditioner.hpp life/lifealg/algebraic_facto.hpp life/lifealg/SolverBase.hpp life/lifealg/triDiagLU.hpp life/lifealg/SolverPETSC.hpp life/lifealg/clapack.hpp life/lifealg/preconditioner.hpp life/lifealg/triDiagCholesky.hpp life/lifealg/dataNewton.hpp life/lifealg/PressureMatrixSolver.hpp life/lifealg/SolverTrilinos.hpp life/lifealg/linesearch_cubic.hpp life/lifealg/generalizedAitken.hpp life/lifealg/inexactLU.hpp life/lifealg/SolverUMFPACK.hpp life/lifealg/cblas.hpp life/lifealg/dataAztec.hpp life/lifealg/EpetraMap.hpp life/lifealg/nonLinRichardson.hpp life/lifealg/linesearch_parabolic.hpp
LIFEMESH_HDRS=
LIFEMESH_HDRS=life/lifemesh/identifier.hpp life/lifemesh/basisElSh.hpp life/lifemesh/geo0D.hpp life/lifemesh/bareItems.hpp life/lifemesh/dataMesh.hpp life/lifemesh/markers.hpp life/lifemesh/partitionMesh.hpp life/lifemesh/basicOneDMesh.hpp life/lifemesh/geoElement.hpp life/lifemesh/regionMesh1D.hpp life/lifemesh/markers_base.hpp life/lifemesh/subDomainNeighbors.hpp life/lifemesh/regionMesh3D.hpp life/lifemesh/geoND.hpp life/lifemesh/regionMesh2D.hpp life/lifemesh/meshEntity.hpp life/lifemesh/mesh_util_base.hpp
LIFEFEM_HDRS=
LIFEFEM_HDRS=life/lifefem/geoMapDG.hpp life/lifefem/oneDBCFunctions.hpp life/lifefem/currentIFDG.hpp life/lifefem/currentHdivFE.hpp life/lifefem/subelements.hpp life/lifefem/elemOper2Fluids.hpp life/lifefem/bdfNS.hpp life/lifefem/refFE.hpp life/lifefem/refEle.hpp life/lifefem/bcFunction.hpp life/lifefem/dofDG.hpp life/lifefem/FESpace.hpp life/lifefem/bcCond.hpp life/lifefem/dofInterface3Dto2D.hpp life/lifefem/dof.hpp life/lifefem/elemOper_ext.hpp life/lifefem/assembGeneric.hpp life/lifefem/currentFE.hpp life/lifefem/staticBdFE.hpp life/lifefem/currentFEDG.hpp life/lifefem/dofInterface3Dto3D.hpp life/lifefem/dataTime.hpp life/lifefem/quadRule.hpp life/lifefem/geoMap.hpp life/lifefem/bdf.hpp life/lifefem/currentBFDG.hpp life/lifefem/oneDBCHandler.hpp life/lifefem/bcHandler.hpp life/lifefem/dataTransient.hpp life/lifefem/dofInterfaceHandler.hpp life/lifefem/Operator01.hpp life/lifefem/assemb.hpp life/lifefem/bcVector.hpp life/lifefem/localDofPattern.hpp life/lifefem/currentBdDG.hpp life/lifefem/refHdivFE.hpp life/lifefem/bdf_template.hpp life/lifefem/dofOneD.hpp life/lifefem/refFEDG.hpp life/lifefem/postProc.hpp life/lifefem/sobolevNorms.hpp life/lifefem/bcManage.hpp life/lifefem/dofByFace.hpp life/lifefem/elemOper.hpp life/lifefem/interpolate.hpp life/lifefem/bdfNS_template.hpp life/lifefem/refHybridFE.hpp life/lifefem/refEleDG.hpp life/lifefem/v2elemVec.hpp life/lifefem/dofInterfaceBase.hpp life/lifefem/regionMesh3D_ALE.hpp life/lifefem/currentBdFE.hpp life/lifefem/values.hpp
LIFESOLVER_HDRS=
LIFESOLVER_HDRS=life/lifesolver/convDiffReactHandler.hpp life/lifesolver/NavierStokesAleSolverCT.hpp life/lifesolver/nsipterms.hpp life/lifesolver/monodomainSolver.hpp life/lifesolver/AFSolvers.hpp life/lifesolver/FSISolver.hpp life/lifesolver/dataBidomain.hpp life/lifesolver/FSIOperator.hpp life/lifesolver/dataIonic.hpp life/lifesolver/oneDModelSolver.hpp life/lifesolver/darcySolverBase.hpp life/lifesolver/ipStabilization.hpp life/lifesolver/dataConvDiffReact.hpp life/lifesolver/oneDNet.hpp life/lifesolver/NavierStokesWithFlux.hpp life/lifesolver/nsip.hpp life/lifesolver/NavierStokesAleSolverPC.hpp life/lifesolver/HarmonicExtensionSolver.hpp life/lifesolver/oneDNonLinModelParam.hpp life/lifesolver/dataNS2Fluids.hpp life/lifesolver/fixedPointBase.hpp life/lifesolver/bidomainSolver.hpp life/lifesolver/NavierStokesHandler.hpp life/lifesolver/dataSimplifiedStructure.hpp life/lifesolver/dataMonodomain.hpp life/lifesolver/oneDModelHandler.hpp life/lifesolver/ionicSolver.hpp life/lifesolver/darcyHandler.hpp life/lifesolver/HyperbolicSolverIP.hpp life/lifesolver/dataDarcy.hpp life/lifesolver/NavierStokesAleSolver.hpp life/lifesolver/LevelSetSolver.hpp life/lifesolver/OseenShapeDerivative.hpp life/lifesolver/darcySolver.hpp life/lifesolver/LevelSetSolverUtils.hpp life/lifesolver/NavierStokesAleHandler.hpp life/lifesolver/exactJacobianBase.hpp life/lifesolver/NavierStokesSolverPC.hpp life/lifesolver/dataElasticStructure.hpp life/lifesolver/simplifiedStructure.hpp life/lifesolver/dataNavierStokes.hpp life/lifesolver/sdStabilization.hpp life/lifesolver/steklovPoincareBase.hpp life/lifesolver/Oseen.hpp life/lifesolver/dataOneDModel.hpp life/lifesolver/VenantKirchhofSolver.hpp life/lifesolver/convDiffReactSolverPC.hpp life/lifesolver/timeSolver.hpp life/lifesolver/vectorFunction1D.hpp life/lifesolver/fluidToMaster.hpp life/lifesolver/parabolicSolver.hpp life/lifesolver/NavierStokesSolver.hpp life/lifesolver/ElasticStructureHandler.hpp life/lifesolver/reducedLinFluid.hpp life/lifesolver/NSSolver2FluidsMixed.hpp
LIFEFILTERS_HDRS=
LIFEFILTERS_HDRS=life/lifefilters/exporter.hpp life/lifefilters/hdf5exporter.hpp life/lifefilters/medit.hpp life/lifefilters/postProc.hpp life/lifefilters/gracePlot.hpp life/lifefilters/gmv_wrtrs.hpp life/lifefilters/readMesh2D.hpp life/lifefilters/openDX_wrtrs.hpp life/lifefilters/selectMarker.hpp life/lifefilters/importer.hpp life/lifefilters/medit_wrtrs.hpp life/lifefilters/ensight7Writer.hpp life/lifefilters/mesh_util.hpp life/lifefilters/gmsh_wrtrs.hpp life/lifefilters/vtk_wrtrs.hpp life/lifefilters/ensight.hpp life/lifefilters/readMesh3D.hpp
nobase_include_HEADERS=lifeconfig.h\
$(LIFECORE_HDRS)\
......
......@@ -39,6 +39,32 @@ EpetraMap::EpetraMap(int NumGlobalElements,
Comm);
}
//! construct a map with entries [1:lagrangeMultipliers] distributed on all the processors
EpetraMap::EpetraMap(std::vector<int> const& lagrangeMultipliers,
const Epetra_Comm& Comm):
M_repeatedEpetra_Map(),
M_uniqueEpetraMap(),
M_exporter(),
M_importer()
{
int NumGlobalElements(lagrangeMultipliers.size());
int NumMyElements (NumGlobalElements);
std::vector<int> MyGlobalElements(lagrangeMultipliers);
int IndexBase = 1;
for (int i(0); i < NumGlobalElements; ++i)
MyGlobalElements[i] = i;
createMap( NumGlobalElements,
NumMyElements,
&MyGlobalElements[0],
IndexBase,
Comm);
}
EpetraMap::map_ptrtype const &
EpetraMap::getMap( EpetraMapType maptype) const
{
......@@ -64,18 +90,25 @@ EpetraMap::getRootMap( int root) const
the maximum id to consider (in the new count)
eg: offset = 2, maxid = 6;
_epetraMap = [ 0 2 5 7 8 10 1]
this = [ 0 3 5 ]
this = [ 0 3 5 7 ]
if needed, indexBase may be changed (default values < 0 means "same as original map")
*/
EpetraMap::EpetraMap(const Epetra_BlockMap& _blockMap, const int offset, const int maxid )
EpetraMap::EpetraMap(const Epetra_BlockMap& _blockMap, const int offset, const int maxid,
int indexbase)
:
M_repeatedEpetra_Map(),
M_uniqueEpetraMap(),
M_exporter(),
M_importer()
{
if (indexbase < 0) indexbase = _blockMap.IndexBase();
std::vector<int> MyGlobalElements;
int* sourceGlobalElements(_blockMap.MyGlobalElements());
int const maxIdOrig(maxid + offset);
int const startIdOrig(offset + indexbase );
int const endIdOrig (startIdOrig + maxid );
const int maxMyElements = std::min(maxid, _blockMap.NumMyElements());
MyGlobalElements.reserve(maxMyElements);
......@@ -83,13 +116,13 @@ EpetraMap::EpetraMap(const Epetra_BlockMap& _blockMap, const int offset, const i
// We consider that the source Map may not be ordered
for (int i(0); i < _blockMap.NumMyElements(); ++i)
if (sourceGlobalElements[i] < maxIdOrig && sourceGlobalElements[i] >= offset)
if (sourceGlobalElements[i] < endIdOrig && sourceGlobalElements[i] >= startIdOrig)
MyGlobalElements.push_back(sourceGlobalElements[i] - offset);
createMap( -1,
MyGlobalElements.size(),
&MyGlobalElements.front(),
_blockMap.IndexBase(),
indexbase,
_blockMap.Comm() );
......@@ -133,7 +166,8 @@ EpetraMap::operator += (const EpetraMap& _epetraMap)
map.push_back(*pointer);
}
int numGlobalElements = getUniqueMap()->NumGlobalElements();
int numGlobalElements = getUniqueMap()->NumGlobalElements()
+ getUniqueMap()->IndexBase() - _epetraMap.getUniqueMap()->IndexBase();
// std::cout << "NumGlobalElements = " << numGlobalElements << std::endl;
......@@ -170,6 +204,21 @@ EpetraMap::operator += (const EpetraMap& _epetraMap)
return *this;
}
EpetraMap &
EpetraMap::operator += (std::vector<int> const& lagrangeMultipliers)
{
if ( lagrangeMultipliers.size() <= 0)
return *this;
ASSERT(this->getUniqueMap(), "operator+=(const int) works only for an existing EpetraMap");
EpetraMap lagrMap(lagrangeMultipliers, Comm());
this->operator+=(lagrMap);
return *this;
}
EpetraMap::EpetraMap(const EpetraMap& _epetraMap):
M_repeatedEpetra_Map(),
......@@ -188,11 +237,19 @@ EpetraMap::createMap(int NumGlobalElements,
int IndexBase,
const Epetra_Comm &Comm)
{
M_repeatedEpetra_Map.reset( new Epetra_Map(NumGlobalElements,
NumMyElements,
MyGlobalElements,
IndexBase,
Comm) );
if (NumMyElements !=0 && MyGlobalElements == 0) // linearMap
M_repeatedEpetra_Map.reset( new Epetra_Map(NumGlobalElements,
NumMyElements,
IndexBase,
Comm) );
else // classic LifeV map
M_repeatedEpetra_Map.reset( new Epetra_Map(NumGlobalElements,
NumMyElements,
MyGlobalElements,
IndexBase,
Comm) );
uniqueMap();
}
......@@ -349,6 +406,5 @@ EpetraMap::MapsAreSimilar( EpetraMap const& _epetraMap) const
}
} // end namespace LifeV
......@@ -70,12 +70,20 @@ public:
//@}
EpetraMap();
// epetra map constructor. To define a linear map, set MyGlobalElements = 0
EpetraMap(int NumGlobalElements,
int NumMyElements,
int* MyGlobalElements,
int IndexBase,
const Epetra_Comm& Comm);
//! construct a map with entries lagrangeMultipliers.
//! Repeated lagrange multipliers will be repeated on the repeatedmap
//! Again: it is not necessary that the lagrangeMltiplier vector is the same on all
//! processors nor that it is different
EpetraMap(std::vector<int> const& lagrangeMultipliers,
const Epetra_Comm& Comm);
// Calls createImportExport from setUp()
template<typename Mesh>
EpetraMap(const RefFE& refFE,
......@@ -92,8 +100,14 @@ public:
/*! Builds a submap of map _epetraMap with a given positive offset and
the maximum id to consider
eg: offset = 2, maxid = 6;
_epetraMap = [ 0 2 5 7 8 10 1]
this = [ 0 3 5 7 ]
if needed, indexBase may be changed (default values < 0 means "same as original map")
*/
EpetraMap(const Epetra_BlockMap& _blockMap, const int offset, const int maxid );
EpetraMap(const Epetra_BlockMap& _blockMap, const int offset, const int maxid,
int indexbase = -1);
~EpetraMap() {}
......@@ -109,6 +123,16 @@ public:
return map;
}
EpetraMap& operator += (std::vector<int> const& lagrangeMultipliers);
EpetraMap operator + (std::vector<int> const& lagrangeMultipliers)
{
EpetraMap map( *this );
map += lagrangeMultipliers;
createImportExport();
return map;
}
Epetra_Comm const& Comm() const { return M_uniqueEpetraMap->Comm(); }
// Epetra_Map* getRepeatedEpetra_Map(){return M_repeatedEpetra_Map;}
......
......@@ -78,6 +78,7 @@ int nonLinRichardson( VectorType& sol,
//----------------------------------------------------------------------
bool const verbose(sol.Comm().MyPID() == 0);
int iter = 0;
......@@ -88,10 +89,12 @@ int nonLinRichardson( VectorType& sol,
Real normResOld = 1;
std::cout << "------------------------------------------------------------------" << std::endl;
std::cout << " NonLinRichardson: starting " << std::endl;
std::cout << "------------------------------------------------------------------" << std::endl;
if (verbose)
{
std::cout << "------------------------------------------------------------------" << std::endl;
std::cout << " NonLinRichardson: starting " << std::endl;
std::cout << "------------------------------------------------------------------" << std::endl;
}
f.evalResidual( residual, sol, iter );
Real normRes = residual.NormInf();
......@@ -116,13 +119,16 @@ int nonLinRichardson( VectorType& sol,
while ( normRes > stop_tol && iter < maxit )
{
std::cout << std::endl;
std::cout << "------------------------------------------------------------------" << std::endl;
std::cout << " NonLinRichardson: iter = " << iter
<< ", residual = " << normRes
<< ", stoping tolerance = " << stop_tol << std::endl;
std::cout << "------------------------------------------------------------------" << std::endl;
std::cout << std::endl;
if (verbose)
{
std::cout << std::endl;
std::cout << "------------------------------------------------------------------" << std::endl;
std::cout << " NonLinRichardson: iter = " << iter
<< ", residual = " << normRes
<< ", stoping tolerance = " << stop_tol << std::endl;
std::cout << "------------------------------------------------------------------" << std::endl;
std::cout << std::endl;
}
iter++;
......@@ -178,24 +184,28 @@ int nonLinRichardson( VectorType& sol,
linearRelTol = std::min<Real>( eta_max,
std::max<Real>( linearRelTol,
.5 * stop_tol / normRes ) );
std::cout << " Newton: forcing term eta = " << linearRelTol << std::endl;
if (verbose)
std::cout << " Newton: forcing term eta = " << linearRelTol << std::endl;
}
}
if ( normRes > stop_tol )
{
std::cout << "!!! NonLinRichardson: convergence fails" << std::endl;
if (verbose)
std::cout << "!!! NonLinRichardson: convergence fails" << std::endl;
maxit = iter;
return 1;
}
//f.displacementOnInterface();
std::cout << "------------------------------------------------------------------" << std::endl;
std::cout << "--- NonLinRichardson: convergence (" << normRes
<<") in " << iter << " iterations\n\n";
std::cout << "------------------------------------------------------------------" << std::endl;
if (verbose)
{
std::cout << "------------------------------------------------------------------" << std::endl;
std::cout << "--- NonLinRichardson: convergence (" << normRes
<<") in " << iter << " iterations\n\n";
std::cout << "------------------------------------------------------------------" << std::endl;
}
maxit = iter;
return 0;
......
......@@ -94,10 +94,10 @@ public:
// Epetra_FECrsMatrix const & getEpetraMatrix(){return M_epetraCrs;}
//! set entries (rVec(i),rVec(i)) to coeff and rest of row r(i) to zero
void diagonalize ( std::vector<UInt> rVec, DataType const coeff );
void diagonalize ( std::vector<UInt> rVec, DataType const coeff, UInt offset=0 );
//! set entry (r,r) to coeff and rest of row r to zero
void diagonalize ( UInt const r, DataType const coeff );
void diagonalize ( UInt const r, DataType const coeff, UInt offset=0 );
/*! apply constraint on all rows rVec
* @param rVec vector of rows
......@@ -106,9 +106,10 @@ public:
* @param datumVec vector of values to constrain entry r of the solution at
*/
void diagonalize( std::vector<UInt> rVec,
DataType const coeff,
vector_type &b,
std::vector<DataType> datumVec );
DataType const coeff,
vector_type &b,
std::vector<DataType> datumVec,
UInt offset=0 );
/*! apply constraint on row r
* @param r row number
* @param coeff value to set entry (r,r) at
......@@ -116,7 +117,8 @@ public:
* @param datum value to constrain entry r of the solution at
*/
void diagonalize( UInt const r, DataType const coeff, vector_type &b,
DataType datum );
DataType datum,
UInt offset=0 );
int getMeanNumEntries() const ;
......@@ -322,7 +324,8 @@ int EpetraMatrix<DataType>::GlobalAssemble()
{
if ( M_epetraCrs.Filled ())
{
std::cout << "Matrix is already filled" << std::endl;
if (M_epetraCrs.Comm().MyPID() == 0)
std::cout << "Matrix is already filled" << std::endl;
return -1;
}
......@@ -338,7 +341,8 @@ void EpetraMatrix<DataType>::insertZeroDiagonal()
if ( M_epetraCrs.Filled ())
{
std::cout << "Matrix is already filled, it is impossible to insert the diagonal now" << std::endl;
if (M_epetraCrs.Comm().MyPID() == 0)
std::cout << "Matrix is already filled, it is impossible to insert the diagonal now" << std::endl;
return;
}
......@@ -359,7 +363,8 @@ void EpetraMatrix<DataType>::insertZeroDiagonal()
//! set entries (rVec(i),rVec(i)) to coeff and rest of row r(i) to zero
template <typename DataType>
void EpetraMatrix<DataType>::diagonalize ( std::vector<UInt> rVec,
DataType const coeff )
DataType const coeff,
UInt offset)
{
const Epetra_Comm& Comm(M_epetraCrs.Comm());
......@@ -395,7 +400,7 @@ void EpetraMatrix<DataType>::diagonalize ( std::vector<UInt> rVec,
Comm.Broadcast(r, sizeVec, p);
for (i=0; i < sizeVec; i++)
diagonalize( Ur[i], coeff);
diagonalize( Ur[i], coeff, offset);
if ( p != MyPID )
{
......@@ -409,7 +414,8 @@ void EpetraMatrix<DataType>::diagonalize ( std::vector<UInt> rVec,
//! set entry (r,r) to coeff and rest of row r to zero
template <typename DataType>
void EpetraMatrix<DataType>::diagonalize( UInt const r,
DataType const coeff)
DataType const coeff,
UInt offset)
{
if ( !M_epetraCrs.Filled() )
......@@ -421,10 +427,10 @@ void EpetraMatrix<DataType>::diagonalize( UInt const r,
const Epetra_Map& colMap(M_epetraCrs.ColMap());
int myCol = colMap.LID(r + 1);
int myCol = colMap.LID(r + 1 + offset);
// row: if r is mine, zero out values
int myRow = rowMap.LID(r + 1);
int myRow = rowMap.LID(r + 1 + offset);
if (myRow >= 0) // I have this row
{
......@@ -459,7 +465,8 @@ template <typename DataType>
void EpetraMatrix<DataType>::diagonalize( std::vector<UInt> rVec,
DataType const coeff,
vector_type &b,
std::vector<DataType> datumVec )
std::vector<DataType> datumVec,
UInt offset)
{
const Epetra_Comm& Comm(M_epetraCrs.Comm());
......@@ -503,7 +510,7 @@ void EpetraMatrix<DataType>::diagonalize( std::vector<UInt> rVec,
Comm.Broadcast(datum,sizeVec, p);
for (i=0; i < sizeVec; i++)
diagonalize( Ur[i], coeff, b, datum[i]);
diagonalize( Ur[i], coeff, b, datum[i], offset);
if ( p != MyPID )
{
......@@ -527,7 +534,8 @@ template <typename DataType>
void EpetraMatrix<DataType>::diagonalize( UInt const r,
DataType const coeff,
vector_type &b,
DataType datum )
DataType datum,
UInt offset)
{
if ( !M_epetraCrs.Filled() )
......@@ -539,7 +547,7 @@ void EpetraMatrix<DataType>::diagonalize( UInt const r,
const Epetra_Map& colMap(M_epetraCrs.ColMap());
int myCol = colMap.LID(r + 1);
int myCol = colMap.LID(r + 1 + offset);
#ifdef EPETRAMATRIX_SYMMETRIC_DIAGONALIZE
if (myCol >= 0) // I have this column
......@@ -554,7 +562,7 @@ void EpetraMatrix<DataType>::diagonalize( UInt const r,
#endif
// row: if r is mine, zero out values
int myRow = rowMap.LID(r + 1);
int myRow = rowMap.LID(r + 1 + offset);
if (myRow >= 0) // I have this row
{
......@@ -572,7 +580,7 @@ void EpetraMatrix<DataType>::diagonalize( UInt const r,
DataType coeff_(coeff);
M_epetraCrs.ReplaceMyValues(myRow, 1, &coeff_, &myCol); // A(r,r) = coeff
b[ r + 1 ] = coeff * datum; // correct right hand side for row r // BASEINDEX + 1
b[ r + 1 + offset] = coeff * datum; // correct right hand side for row r // BASEINDEX + 1
}
......
......@@ -50,14 +50,14 @@ EpetraVector::EpetraVector( const EpetraVector& _vector):
EpetraVector::EpetraVector( const EpetraMap& _map, EpetraMapType maptype ):
M_epetraVector(*_map.getMap(maptype), false),
M_epetraVector(*_map.getMap(maptype)),
M_epetraMap (new EpetraMap(_map)),
M_maptype (maptype)
{
}
EpetraVector::EpetraVector( const EpetraVector& _vector, EpetraMapType maptype):
M_epetraVector(*_vector.M_epetraMap->getMap(maptype), false),
M_epetraVector(*_vector.M_epetraMap->getMap(maptype)),
M_epetraMap (_vector.M_epetraMap),
M_maptype (maptype)
{
......@@ -66,7 +66,7 @@ EpetraVector::EpetraVector( const EpetraVector& _vector, EpetraMapType maptype):
EpetraVector::EpetraVector( const EpetraVector& _vector, EpetraMapType maptype,
Epetra_CombineMode combineMode):
M_epetraVector(*_vector.M_epetraMap->getMap(maptype), false),
M_epetraVector(*_vector.M_epetraMap->getMap(maptype)),
M_epetraMap (_vector.M_epetraMap),
M_maptype (maptype)
{
......@@ -102,7 +102,7 @@ EpetraVector::EpetraVector( const Epetra_MultiVector& _vector,
boost::shared_ptr<EpetraMap> _map,
EpetraMapType maptype )
:
M_epetraVector( *_map->getMap(maptype), false),
M_epetraVector( *_map->getMap(maptype)),
M_epetraMap ( _map),
M_maptype (maptype)
{
......@@ -116,7 +116,7 @@ EpetraVector::EpetraVector( const Epetra_MultiVector& _vector,
// Copies _vector to a vector which resides only on the processor "reduceToProc"
EpetraVector::EpetraVector( const EpetraVector& _vector, const int reduceToProc):
M_epetraVector(_vector.M_epetraMap->getRootMap(reduceToProc), false),
M_epetraVector(_vector.M_epetraMap->getRootMap(reduceToProc)),
M_epetraMap (),
M_maptype (Unique)
{
......@@ -198,9 +198,9 @@ int EpetraVector::checkLID(const UInt row) const
//! if row is mine sets this[row] = value and return true
//! if row is not mine and if the numCpus > 1, returns false
//! if row is not mine and if the numCpus == 1, asserts
bool EpetraVector::checkAndSet(const UInt row, const data_type& value)
bool EpetraVector::checkAndSet(const UInt row, const data_type& value, UInt offset)
{
int lrow = checkLID(row);
int lrow = checkLID(row + offset);
if (lrow < 0)
return false;
......@@ -409,14 +409,12 @@ operator = (const EpetraVector& _vector)
return *this;
}
}
// We hope we are guessing right
switch (M_maptype) {
case Unique:
//
return Export(_vector.M_epetraVector, Add);
switch (_vector.M_maptype) {
case Repeated:
if (_vector.M_maptype != Repeated)
//
if (M_maptype != Repeated)
return Export(_vector.M_epetraVector, Add);
case Unique:
return Import(_vector.M_epetraVector, Add);
}
......
......@@ -110,7 +110,7 @@ public:
//! if row is mine sets this[row] = value and return true
//! if row is not mine and if the numCpus > 1, returns false
//! if row is not mine and if the numCpus == 1, asserts
bool checkAndSet(const UInt row, const data_type& value);
bool checkAndSet(const UInt row, const data_type& value, UInt offset=0);
//! Set the row row of the vector to value. If it isn't on this processor,
//! store it and send it and send it at next GlobalAssemble
......
......@@ -31,7 +31,7 @@ namespace LifeV
/#Usage:
/ SimpleVect<float> a(10);
/ a(10)=90; // a[9] will contain 90.0
/ SimpleArray,int> b(3,5) // an arrray with 3 rows
/ SimpleArray<int> b(3,5) // an arrray with 3 rows
/ and 5 columns
/ b(3,2)=5;
/ b.reshape(2,3) // now b is 2x3
......@@ -206,7 +206,7 @@ template <typename T, int OFFSETVEC>
SimpleArray<T, OFFSETVEC>::SimpleArray( size_type ntot )
:
std::vector<T>( ntot ),
_M_nrows( nrows ),
_M_nrows( ntot ),
_M_ncols( 1 )
{}
......
......@@ -69,6 +69,11 @@ public:
_mat = 0.0;
};
void showMe( std::ostream& c = std::cout );
void operator *= (Real coef)
{
this->_mat *= coef;
}
private:
matrix_type _mat;
......
......@@ -39,7 +39,7 @@
#include <life/lifefem/refFE.hpp>
#include <life/lifefem/dof.hpp>
#include <life/lifefem/geoMap.hpp>
//#include <life/lifefem/bcHandler.hpp>
#include <life/lifefem/bcHandler.hpp>
#include <life/lifearray/pattern.hpp>
#include <cmath>
#include <sstream>
......@@ -50,6 +50,7 @@
#include <life/lifemesh/partitionMesh.hpp>
#include <life/lifefem/currentFE.hpp>
#include <life/lifefem/currentBdFE.hpp>
#include <life/lifefem/sobolevNorms.hpp>
using std::pair;
......@@ -180,8 +181,10 @@ public:
template<typename vector_type>
void interpolate( const Function& fct, vector_type& vect, const Real time = 0. );
//! calculate L2 pressure error for given exact pressure function
//! takes into account a possible offset by a constant
template<typename vector_type>