From: vbd Date: Wed, 13 Jul 2011 08:13:11 +0000 (+0000) Subject: sources au 13/07 - ok cells pb faces X-Git-Tag: MEDPartitioner_first_compilable_version~11 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=6c60e07d2415d2ae4397cb18a236b28b4589a62f;p=tools%2Fmedcoupling.git sources au 13/07 - ok cells pb faces --- diff --git a/src/INTERP_KERNEL/Bases/NormalizedUnstructuredMesh.hxx b/src/INTERP_KERNEL/Bases/NormalizedUnstructuredMesh.hxx index 5d2350cae..5ccfdbf38 100644 --- a/src/INTERP_KERNEL/Bases/NormalizedUnstructuredMesh.hxx +++ b/src/INTERP_KERNEL/Bases/NormalizedUnstructuredMesh.hxx @@ -31,7 +31,7 @@ namespace INTERP_KERNEL typedef enum { - NORM_POINT0 = 0, + NORM_POINT1 = 0, NORM_SEG2 = 1, NORM_SEG3 = 2, NORM_TRI3 = 3, diff --git a/src/INTERP_KERNEL/CellModel.cxx b/src/INTERP_KERNEL/CellModel.cxx index add613149..747419dff 100644 --- a/src/INTERP_KERNEL/CellModel.cxx +++ b/src/INTERP_KERNEL/CellModel.cxx @@ -27,7 +27,7 @@ namespace INTERP_KERNEL { - const char *CellModel::CELL_TYPES_REPR[]={"NORM_POINT0", "NORM_SEG2", "NORM_SEG3", "NORM_TRI3", "NORM_QUAD4",// 0->4 + const char *CellModel::CELL_TYPES_REPR[]={"NORM_POINT1", "NORM_SEG2", "NORM_SEG3", "NORM_TRI3", "NORM_QUAD4",// 0->4 "NORM_POLYGON", "NORM_TRI6", "" , "NORM_QUAD8", "",//5->9 "", "", "", "", "NORM_TETRA4",//10->14 "NORM_PYRA5", "NORM_PENTA6", "", "NORM_HEXA8", "",//15->19 @@ -78,7 +78,7 @@ namespace INTERP_KERNEL void CellModel::buildUniqueInstance() { - _map_of_unique_instance.insert(std::make_pair(NORM_POINT0,CellModel(NORM_POINT0))); + _map_of_unique_instance.insert(std::make_pair(NORM_POINT1,CellModel(NORM_POINT1))); _map_of_unique_instance.insert(std::make_pair(NORM_SEG2,CellModel(NORM_SEG2))); _map_of_unique_instance.insert(std::make_pair(NORM_SEG3,CellModel(NORM_SEG3))); _map_of_unique_instance.insert(std::make_pair(NORM_TRI3,CellModel(NORM_TRI3))); @@ -106,9 +106,9 @@ namespace INTERP_KERNEL _quadratic_type=NORM_ERROR; switch(type) { - case NORM_POINT0: + case NORM_POINT1: { - _nb_of_pts=0; _nb_of_sons=0; _dim=0; _extruded_type=NORM_SEG2; _is_simplex=true; + _nb_of_pts=1; _nb_of_sons=0; _dim=0; _extruded_type=NORM_SEG2; _is_simplex=true; } break; case NORM_SEG2: diff --git a/src/INTERP_KERNEL/ExprEval/InterpKernelAsmX86.cxx b/src/INTERP_KERNEL/ExprEval/InterpKernelAsmX86.cxx index ee7116fd2..479b5c78e 100644 --- a/src/INTERP_KERNEL/ExprEval/InterpKernelAsmX86.cxx +++ b/src/INTERP_KERNEL/ExprEval/InterpKernelAsmX86.cxx @@ -23,6 +23,14 @@ #include #include +#ifdef _POSIX_MAPPED_FILES +#include +#else +#ifdef WNT +#include +#endif +#endif + const char *INTERP_KERNEL::AsmX86::OPS[NB_OF_OPS]={"mov","push","pop","fld","faddp","fsubp","fmulp","fdivp","fcos","fsin","fabs","fchs","fsqrt","sub","add","ret","leave","movsd","fst"}; std::vector INTERP_KERNEL::AsmX86::convertIntoMachineLangage(const std::vector& asmb) const throw(INTERP_KERNEL::Exception) @@ -33,11 +41,20 @@ std::vector INTERP_KERNEL::AsmX86::convertIntoMachineLangage(const std::ve return ret; } -char *INTERP_KERNEL::AsmX86::convertMachineLangageInBasic(const std::vector& ml, int& lgth) const +char *INTERP_KERNEL::AsmX86::copyToExecMemZone(const std::vector& ml, unsigned& offset) const { - lgth=ml.size(); - char *ret=new char[lgth]; - std::copy(ml.begin(),ml.end(),ret); + char *ret=0; + int lgth=ml.size(); +#ifdef _POSIX_MAPPED_FILES + ret=(char *)mmap(0,lgth,PROT_EXEC | PROT_WRITE,MAP_ANONYMOUS | MAP_PRIVATE,-1,0); +#else +#ifdef WNT + HANDLE h=CreateFileMapping(INVALID_HANDLE_VALUE,NULL,PAGE_EXECUTE_READWRITE,0,lgth,NULL); + ret=(char *)MapViewOfFile(h,FILE_MAP_EXECUTE | FILE_MAP_READ | FILE_MAP_WRITE,0,0,lgth); +#endif +#endif + if(ret) + std::copy(ml.begin(),ml.end(),ret); return ret; } diff --git a/src/INTERP_KERNEL/ExprEval/InterpKernelAsmX86.hxx b/src/INTERP_KERNEL/ExprEval/InterpKernelAsmX86.hxx index 11b8733e5..fd7bc9afc 100644 --- a/src/INTERP_KERNEL/ExprEval/InterpKernelAsmX86.hxx +++ b/src/INTERP_KERNEL/ExprEval/InterpKernelAsmX86.hxx @@ -32,7 +32,7 @@ namespace INTERP_KERNEL { public: std::vector convertIntoMachineLangage(const std::vector& asmb) const throw(INTERP_KERNEL::Exception); - char *convertMachineLangageInBasic(const std::vector& ml, int& lgth) const; + char *copyToExecMemZone(const std::vector& ml, unsigned& offset) const; private: void convertOneInstructionInML(const std::string& inst, std::vector& ml) const throw(INTERP_KERNEL::Exception); private: diff --git a/src/INTERP_KERNEL/ExprEval/InterpKernelExprParser.cxx b/src/INTERP_KERNEL/ExprEval/InterpKernelExprParser.cxx index 585fac2f8..1f6b497b5 100644 --- a/src/INTERP_KERNEL/ExprEval/InterpKernelExprParser.cxx +++ b/src/INTERP_KERNEL/ExprEval/InterpKernelExprParser.cxx @@ -28,14 +28,6 @@ #include #include -#ifdef _POSIX_MAPPED_FILES -#include -#else -#ifdef WNT -#include -#endif -#endif - using namespace INTERP_KERNEL; const char LeafExprVar::END_OF_RECOGNIZED_VAR[]="Vec"; @@ -719,21 +711,8 @@ char *ExprParser::compileX86() const for(std::vector::const_iterator iter=output.begin();iter!=output.end();iter++) std::cout << std::hex << (int)((unsigned char)(*iter)) << " "; std::cout << std::endl; - int lgth; - char *lm=asmb.convertMachineLangageInBasic(output,lgth); - char *ret=0; -#ifdef _POSIX_MAPPED_FILES - ret=(char *)mmap(0,lgth,PROT_EXEC | PROT_WRITE,MAP_ANONYMOUS | MAP_PRIVATE,-1,0); -#else -#ifdef WNT - HANDLE h=CreateFileMapping(INVALID_HANDLE_VALUE,NULL,PAGE_EXECUTE_READWRITE,0,lgth,NULL); - ret=(char *)MapViewOfFile(h,FILE_MAP_EXECUTE | FILE_MAP_READ | FILE_MAP_WRITE,0,0,lgth); -#endif -#endif - if(ret) - std::copy(lm,lm+lgth,ret); - delete [] lm; - return ret; + unsigned offset; + return asmb.copyToExecMemZone(output,offset); } char *ExprParser::compileX86_64() const @@ -757,21 +736,8 @@ char *ExprParser::compileX86_64() const for(std::vector::const_iterator iter=output.begin();iter!=output.end();iter++) std::cout << std::hex << (int)((unsigned char)(*iter)) << " "; std::cout << std::endl; - int lgth; - char *lm=asmb.convertMachineLangageInBasic(output,lgth); - char *ret=0; -#ifdef _POSIX_MAPPED_FILES - ret=(char *)mmap(0,lgth,PROT_EXEC | PROT_WRITE,MAP_ANONYMOUS | MAP_PRIVATE,-1,0); -#else -#ifdef WNT - HANDLE h=CreateFileMapping(INVALID_HANDLE_VALUE,NULL,PAGE_EXECUTE_READWRITE,0,lgth,NULL); - ret=(char *)MapViewOfFile(h,FILE_MAP_EXECUTE | FILE_MAP_READ | FILE_MAP_WRITE,0,0,lgth); -#endif -#endif - if(ret) - std::copy(lm,lm+lgth,ret); - delete [] lm; - return ret; + unsigned offset; + return asmb.copyToExecMemZone(output,offset); } void ExprParser::compileX86LowLev(std::vector& ass) const @@ -782,9 +748,9 @@ void ExprParser::compileX86LowLev(std::vector& ass) const { for(std::list::const_iterator iter=_sub_expr.begin();iter!=_sub_expr.end();iter++) (*iter).compileX86LowLev(ass); - for(std::list::const_iterator iter2=_func_btw_sub_expr.begin();iter2!=_func_btw_sub_expr.end();iter2++) - (*iter2)->operateX86(ass); } + for(std::list::const_iterator iter2=_func_btw_sub_expr.begin();iter2!=_func_btw_sub_expr.end();iter2++) + (*iter2)->operateX86(ass); } void ExprParser::compileX86_64LowLev(std::vector& ass) const @@ -795,9 +761,9 @@ void ExprParser::compileX86_64LowLev(std::vector& ass) const { for(std::list::const_iterator iter=_sub_expr.begin();iter!=_sub_expr.end();iter++) (*iter).compileX86_64LowLev(ass); - for(std::list::const_iterator iter2=_func_btw_sub_expr.begin();iter2!=_func_btw_sub_expr.end();iter2++) - (*iter2)->operateX86(ass); } + for(std::list::const_iterator iter2=_func_btw_sub_expr.begin();iter2!=_func_btw_sub_expr.end();iter2++) + (*iter2)->operateX86(ass); } void LeafExprVal::compileX86(std::vector& ass) const diff --git a/src/INTERP_KERNEL/GaussPoints/GaussCoords.cxx b/src/INTERP_KERNEL/GaussPoints/GaussCoords.cxx new file mode 100644 index 000000000..73f3b4efc --- /dev/null +++ b/src/INTERP_KERNEL/GaussPoints/GaussCoords.cxx @@ -0,0 +1,1258 @@ +// Copyright (C) 2007-2010 CEA/DEN, EDF R&D +// +// 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. +// +// 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 +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +//Local includes +#include "GaussCoords.hxx" +#include "CellModel.hxx" +using namespace INTERP_KERNEL; + +//STL includes +#include +#include +#include + +//#define MYDEBUG + +#ifdef MYDEBUG +#include +#endif + +//Define common part of the code in the MACRO +//--------------------------------------------------------------- +#define LOCAL_COORD_MACRO_BEGIN \ +myLocalReferenceCoord.resize( myLocalRefDim*myLocalNbRef ); \ +for( int refId = 0; refId < myLocalNbRef; refId++ ) { \ + double* coords = &myLocalReferenceCoord[ refId*myLocalRefDim ]; \ + switch(refId) { + +//--------------------------------------------------------------- +#define LOCAL_COORD_MACRO_END \ + }\ +} + +//--------------------------------------------------------------- +#define SHAPE_FUN_MACRO_BEGIN \ +for( int gaussId = 0 ; gaussId < myNbGauss ; gaussId++ ) { \ + double* funValue = &myFunctionValue[ gaussId * myNbRef ]; \ + const double* gc = &myGaussCoord[ gaussId * GetGaussCoordDim() ]; + +//--------------------------------------------------------------- +#define SHAPE_FUN_MACRO_END \ +} + +#define CHECK_MACRO \ +if( ! aSatify ) { \ + std::ostringstream stream; \ + stream << "Error in the gauss localization for the cell with type "; \ + stream << cellModel.getRepr(); \ + stream << " !!!"; \ + throw INTERP_KERNEL::Exception(stream.str().c_str()); \ +} + + +//--------------------------------------------------------------- +static bool IsEqual(double theLeft, double theRight) { + static double EPS = 1.0E-3; + if(fabs(theLeft) + fabs(theRight) > EPS) + return fabs(theLeft-theRight)/(fabs(theLeft)+fabs(theRight)) < EPS; + return true; +} + + +//////////////////////////////////////////////////////////////////////////////////////////////// +// GAUSS INFO CLASS // +//////////////////////////////////////////////////////////////////////////////////////////////// + +/*! + * Constructor of the GaussInfo + */ +GaussInfo::GaussInfo( NormalizedCellType theGeometry, + const DataVector& theGaussCoord, + int theNbGauss, + const DataVector& theReferenceCoord, + int theNbRef ) : + myGeometry(theGeometry), + myGaussCoord(theGaussCoord), + myNbGauss(theNbGauss), + myReferenceCoord(theReferenceCoord), + myNbRef(theNbRef) { + + //Allocate shape function values + myFunctionValue.resize( myNbGauss * myNbRef ); +} + +/*! + * Destructor + */ +GaussInfo::~GaussInfo() { } + +/*! + * Return dimension of the gauss coordinates + */ +int GaussInfo::GetGaussCoordDim() const { + if( myNbGauss ) { + return myGaussCoord.size()/myNbGauss; + } else { + return 0; + } +} + +/*! + * Return dimension of the reference coordinates + */ +int GaussInfo::GetReferenceCoordDim() const { + if( myNbRef ) { + return myReferenceCoord.size()/myNbRef; + } else { + return 0; + } +} + +/*! + * Return type of the cell. + */ +NormalizedCellType GaussInfo::GetCellType() const { + return myGeometry; +} + +/*! + * Return Nb of the gauss points. + */ +int GaussInfo::GetNbGauss() const { + return myNbGauss; +} + +/*! + * Return Nb of the reference coordinates. + */ +int GaussInfo::GetNbRef() const { + return myNbRef; +} + +/*! + * Check coordinates + */ +bool GaussInfo::isSatisfy() { + + bool anIsSatisfy = ((myLocalNbRef == myNbRef) && (myLocalRefDim == GetReferenceCoordDim())); + //Check coordinates + if(anIsSatisfy) { + for( int refId = 0; refId < myLocalNbRef; refId++ ) { + double* refCoord = &myReferenceCoord[ refId*myLocalRefDim ]; + double* localRefCoord = &myLocalReferenceCoord[ refId*myLocalRefDim ]; + bool anIsEqual = false; + for( int dimId = 0; dimId < myLocalRefDim; dimId++ ) { + anIsEqual = IsEqual( localRefCoord[dimId], refCoord[dimId]); + if(!anIsEqual ) { + return false; + } + } + } + } + return anIsSatisfy; +} + +/*! + * Initialize the internal vectors + */ +void GaussInfo::InitLocalInfo() throw (INTERP_KERNEL::Exception) { + bool aSatify = false; + const CellModel& cellModel=CellModel::getCellModel(myGeometry); + switch( myGeometry ) + { + case NORM_SEG2: + myLocalRefDim = 1; myLocalNbRef = 2; Seg2Init(); + aSatify = isSatisfy(); + CHECK_MACRO; + break; + + case NORM_SEG3: + myLocalRefDim = 1; myLocalNbRef = 3; Seg3Init(); + aSatify = isSatisfy(); + CHECK_MACRO; + break; + + case NORM_TRI3: + myLocalRefDim = 2; myLocalNbRef = 3; Tria3aInit(); + aSatify = isSatisfy(); + + if(!aSatify) { + Tria3bInit(); + aSatify = isSatisfy(); + CHECK_MACRO; + } + break; + + case NORM_TRI6: + myLocalRefDim = 2; myLocalNbRef = 6; Tria6aInit(); + aSatify = isSatisfy(); + if(!aSatify) { + Tria6bInit(); + aSatify = isSatisfy(); + CHECK_MACRO; + } + break; + + case NORM_QUAD4: + myLocalRefDim = 2; myLocalNbRef = 4; Quad4aInit(); + aSatify = isSatisfy(); + + if(!aSatify) { + Quad4bInit(); + aSatify = isSatisfy(); + CHECK_MACRO; + } + break; + + case NORM_QUAD8: + myLocalRefDim = 2; myLocalNbRef = 8; Quad8aInit(); + aSatify = isSatisfy(); + + if(!aSatify) { + Quad8bInit(); + aSatify = isSatisfy(); + CHECK_MACRO; + } + break; + + case NORM_TETRA4: + myLocalRefDim = 3; myLocalNbRef = 4; Tetra4aInit(); + aSatify = isSatisfy(); + + if(!aSatify) { + Tetra4bInit(); + aSatify = isSatisfy(); + CHECK_MACRO; + } + break; + + case NORM_TETRA10: + myLocalRefDim = 3; myLocalNbRef = 10; Tetra10aInit(); + aSatify = isSatisfy(); + + if(!aSatify) { + Tetra10bInit(); + aSatify = isSatisfy(); + CHECK_MACRO; + } + break; + + case NORM_PYRA5: + myLocalRefDim = 3; myLocalNbRef = 5; Pyra5aInit(); + aSatify = isSatisfy(); + + if(!aSatify) { + Pyra5bInit(); + aSatify = isSatisfy(); + CHECK_MACRO; + } + break; + + case NORM_PYRA13: + myLocalRefDim = 3; myLocalNbRef = 13; Pyra13aInit(); + aSatify = isSatisfy(); + + if(!aSatify) { + Pyra13bInit(); + aSatify = isSatisfy(); + CHECK_MACRO; + } + break; + + case NORM_PENTA6: + myLocalRefDim = 3; myLocalNbRef = 6; Penta6aInit(); + aSatify = isSatisfy(); + + if(!aSatify) { + Penta6bInit(); + aSatify = isSatisfy(); + CHECK_MACRO; + } + break; + + case NORM_PENTA15: + myLocalRefDim = 3; myLocalNbRef = 15; Penta15aInit(); + aSatify = isSatisfy(); + + if(!aSatify) { + Penta15bInit(); + aSatify = isSatisfy(); + CHECK_MACRO; + } + break; + + case NORM_HEXA8: + myLocalRefDim = 3; myLocalNbRef = 8; Hexa8aInit(); + aSatify = isSatisfy(); + + if(!aSatify) { + Hexa8aInit(); + aSatify = isSatisfy(); + CHECK_MACRO; + } + break; + + case NORM_HEXA20: + myLocalRefDim = 3; myLocalNbRef = 20; Hexa20aInit(); + aSatify = isSatisfy(); + + if(!aSatify) { + Hexa20aInit(); + aSatify = isSatisfy(); + CHECK_MACRO; + } + break; + + default: + throw INTERP_KERNEL::Exception("Bad Cell type !!!"); + break; + } +} + +/** + * Return shape function value by node id + */ +const double* GaussInfo::GetFunctionValues( const int theGaussId ) const { + return &myFunctionValue[ myNbRef*theGaussId ]; +} + +/*! + * Init Segment 2 Reference coordinates ans Shape function. + */ +void GaussInfo::Seg2Init() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = -1.0; break; + case 1: coords[0] = 1.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = 0.5*(1.0 - gc[0]); + funValue[1] = 0.5*(1.0 + gc[0]); + SHAPE_FUN_MACRO_END +} + +/*! + * Init Segment 3 Reference coordinates ans Shape function. + */ +void GaussInfo::Seg3Init() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = -1.0; break; + case 1: coords[0] = 1.0; break; + case 2: coords[0] = 0.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = 0.5*(1.0 - gc[0])*gc[0]; + funValue[1] = 0.5*(1.0 + gc[0])*gc[0]; + funValue[2] = (1.0 + gc[0])*(1.0 - gc[0]); + SHAPE_FUN_MACRO_END +} + +/*! + * Init Triangle Reference coordinates ans Shape function. + * Case A. + */ +void GaussInfo::Tria3aInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = -1.0; coords[1] = 1.0; break; + case 1: coords[0] = -1.0; coords[1] = -1.0; break; + case 2: coords[0] = 1.0; coords[1] = -1.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = 0.5*(1.0 + gc[1]); + funValue[1] = -0.5*(gc[0] + gc[1]); + funValue[2] = 0.5*(1.0 + gc[0]); + SHAPE_FUN_MACRO_END +} + +/*! + * Init Triangle Reference coordinates ans Shape function. + * Case B. + */ +void GaussInfo::Tria3bInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = 0.0; coords[1] = 0.0; break; + case 1: coords[0] = 1.0; coords[1] = 0.0; break; + case 2: coords[0] = 0.0; coords[1] = 1.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = 1.0 - gc[0] - gc[1]; + funValue[1] = gc[0]; + funValue[2] = gc[1]; + SHAPE_FUN_MACRO_END +} + +/*! + * Init Quadratic Triangle Reference coordinates ans Shape function. + * Case A. + */ +void GaussInfo::Tria6aInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = -1.0; coords[1] = 1.0; break; + case 1: coords[0] = -1.0; coords[1] = -1.0; break; + case 2: coords[0] = 1.0; coords[1] = -1.0; break; + case 3: coords[0] = -1.0; coords[1] = 1.0; break; + case 4: coords[0] = 0.0; coords[1] = -1.0; break; + case 5: coords[0] = 0.0; coords[1] = 0.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = 0.5*(1.0 + gc[1])*gc[1]; + funValue[1] = 0.5*(gc[0] + gc[1])*(gc[0] + gc[1] + 1); + funValue[2] = 0.5*(1.0 + gc[0])*gc[0]; + funValue[3] = -1.0*(1.0 + gc[1])*(gc[0] + gc[1]); + funValue[4] = -1.0*(1.0 + gc[0])*(gc[0] + gc[1]); + funValue[5] = (1.0 + gc[1])*(1.0 + gc[1]); + SHAPE_FUN_MACRO_END +} + +/*! + * Init Quadratic Triangle Reference coordinates ans Shape function. + * Case B. + */ +void GaussInfo::Tria6bInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = 0.0; coords[1] = 0.0; break; + case 1: coords[0] = 1.0; coords[1] = 0.0; break; + case 2: coords[0] = 0.0; coords[1] = 1.0; break; + case 3: coords[0] = 0.5; coords[1] = 0.0; break; + case 4: coords[0] = 0.5; coords[1] = 0.5; break; + case 5: coords[0] = 0.0; coords[1] = 0.5; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = (1.0 - gc[0] - gc[1])*(1.0 - 2.0*gc[0] - 2.0*gc[1]); + funValue[1] = gc[0]*(2.0*gc[0] - 1.0); + funValue[2] = gc[1]*(2.0*gc[1] - 1.0); + funValue[3] = 4.0*gc[0]*(1.0 - gc[0] - gc[1]); + funValue[4] = 4.0*gc[0]*gc[1]; + funValue[5] = 4.0*gc[1]*(1.0 - gc[0] - gc[1]); + SHAPE_FUN_MACRO_END +} + +/*! + * Init Quadrangle Reference coordinates ans Shape function. + * Case A. + */ +void GaussInfo::Quad4aInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = -1.0; coords[1] = 1.0; break; + case 1: coords[0] = -1.0; coords[1] = -1.0; break; + case 2: coords[0] = 1.0; coords[1] = -1.0; break; + case 3: coords[0] = 1.0; coords[1] = 1.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0]); + funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0]); + funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0]); + funValue[3] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]); + SHAPE_FUN_MACRO_END +} + +/*! + * Init Quadrangle Reference coordinates ans Shape function. + * Case B. + */ +void GaussInfo::Quad4bInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = -1.0; coords[1] = -1.0; break; + case 1: coords[0] = 1.0; coords[1] = -1.0; break; + case 2: coords[0] = 1.0; coords[1] = 1.0; break; + case 3: coords[0] = -1.0; coords[1] = 1.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]); + funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]); + funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]); + funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]); + SHAPE_FUN_MACRO_END +} + + +/*! + * Init Quadratic Quadrangle Reference coordinates ans Shape function. + * Case A. + */ +void GaussInfo::Quad8aInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = -1.0; coords[1] = 1.0; break; + case 1: coords[0] = -1.0; coords[1] = -1.0; break; + case 2: coords[0] = 1.0; coords[1] = -1.0; break; + case 3: coords[0] = 1.0; coords[1] = 1.0; break; + case 4: coords[0] = -1.0; coords[1] = 0.0; break; + case 5: coords[0] = 0.0; coords[1] = -1.0; break; + case 6: coords[0] = 1.0; coords[1] = 0.0; break; + case 7: coords[0] = 0.0; coords[1] = 1.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0])*(gc[1] - gc[0] - 1.0); + funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0])*(-gc[1] - gc[0] - 1.0); + funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0])*(-gc[1] + gc[0] - 1.0); + funValue[3] = 0.25*(1.0 + gc[1])*(1.0 + gc[0])*(gc[1] + gc[0] - 1.0); + funValue[4] = 0.5*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[1]); + funValue[5] = 0.5*(1.0 - gc[1])*(1.0 - gc[0])*(1.0 + gc[0]); + funValue[6] = 0.5*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[1]); + funValue[7] = 0.5*(1.0 + gc[1])*(1.0 - gc[0])*(1.0 + gc[0]); + SHAPE_FUN_MACRO_END +} + +/*! + * Init Quadratic Quadrangle Reference coordinates ans Shape function. + * Case B. + */ +void GaussInfo::Quad8bInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = -1.0; coords[1] = -1.0; break; + case 1: coords[0] = 1.0; coords[1] = -1.0; break; + case 2: coords[0] = 1.0; coords[1] = 1.0; break; + case 3: coords[0] = -1.0; coords[1] = 1.0; break; + case 4: coords[0] = 0.0; coords[1] = -1.0; break; + case 5: coords[0] = 1.0; coords[1] = 0.0; break; + case 6: coords[0] = 0.0; coords[1] = 1.0; break; + case 7: coords[0] = -1.0; coords[1] = 0.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1])*(-1.0 - gc[0] - gc[1]); + funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1])*(-1.0 + gc[0] - gc[1]); + funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1])*(-1.0 + gc[0] + gc[1]); + funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1])*(-1.0 - gc[0] + gc[1]); + funValue[4] = 0.5*(1.0 - gc[0]*gc[0])*(1.0 - gc[1]); + funValue[5] = 0.5*(1.0 - gc[1]*gc[1])*(1.0 + gc[0]); + funValue[6] = 0.5*(1.0 - gc[0]*gc[0])*(1.0 + gc[1]); + funValue[7] = 0.5*(1.0 - gc[1]*gc[1])*(1.0 - gc[0]); + SHAPE_FUN_MACRO_END +} + +/*! + * Init Tetrahedron Reference coordinates ans Shape function. + * Case A. + */ +void GaussInfo::Tetra4aInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = 0.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 1: coords[0] = 0.0; coords[1] = 0.0; coords[2] = 1.0; break; + case 2: coords[0] = 0.0; coords[1] = 0.0; coords[2] = 0.0; break; + case 3: coords[0] = 1.0; coords[1] = 0.0; coords[2] = 0.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = gc[1]; + funValue[1] = gc[2]; + funValue[2] = 1.0 - gc[0] - gc[1] - gc[2]; + funValue[3] = gc[0]; + SHAPE_FUN_MACRO_END +} + +/*! + * Init Tetrahedron Reference coordinates ans Shape function. + * Case B. + */ +void GaussInfo::Tetra4bInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = 0.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 2: coords[0] = 0.0; coords[1] = 0.0; coords[2] = 1.0; break; + case 1: coords[0] = 0.0; coords[1] = 0.0; coords[2] = 0.0; break; + case 3: coords[0] = 1.0; coords[1] = 0.0; coords[2] = 0.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = gc[1]; + funValue[2] = gc[2]; + funValue[1] = 1.0 - gc[0] - gc[1] - gc[2]; + funValue[3] = gc[0]; + SHAPE_FUN_MACRO_END + +} + +/*! + * Init Quadratic Tetrahedron Reference coordinates ans Shape function. + * Case A. + */ +void GaussInfo::Tetra10aInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = 0.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 1: coords[0] = 0.0; coords[1] = 0.0; coords[2] = 1.0; break; + case 2: coords[0] = 0.0; coords[1] = 0.0; coords[2] = 0.0; break; + case 3: coords[0] = 1.0; coords[1] = 0.0; coords[2] = 0.0; break; + case 4: coords[0] = 0.0; coords[1] = 0.5; coords[2] = 0.5; break; + case 5: coords[0] = 0.0; coords[1] = 0.0; coords[2] = 0.5; break; + case 6: coords[0] = 0.0; coords[1] = 0.5; coords[2] = 0.0; break; + case 7: coords[0] = 0.5; coords[1] = 0.5; coords[2] = 0.0; break; + case 8: coords[0] = 0.5; coords[1] = 0.0; coords[2] = 0.5; break; + case 9: coords[0] = 0.5; coords[1] = 0.0; coords[2] = 0.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = gc[1]*(2.0*gc[1] - 1.0); + funValue[1] = gc[2]*(2.0*gc[2] - 1.0); + funValue[2] = (1.0 - gc[0] - gc[1] - gc[2])*(1.0 - 2.0*gc[0] - 2.0*gc[1] - 2.0*gc[2]); + funValue[3] = gc[0]*(2.0*gc[0] - 1.0); + funValue[4] = 4.0*gc[1]*gc[2]; + funValue[5] = 4.0*gc[2]*(1.0 - gc[0] - gc[1] - gc[2]); + funValue[6] = 4.0*gc[1]*(1.0 - gc[0] - gc[1] - gc[2]); + funValue[7] = 4.0*gc[0]*gc[1]; + funValue[8] = 4.0*gc[0]*gc[2]; + funValue[9] = 4.0*gc[0]*(1.0 - gc[0] - gc[1] - gc[2]); + SHAPE_FUN_MACRO_END +} + +/*! + * Init Quadratic Tetrahedron Reference coordinates ans Shape function. + * Case B. + */ +void GaussInfo::Tetra10bInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = 0.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 2: coords[0] = 0.0; coords[1] = 0.0; coords[2] = 1.0; break; + case 1: coords[0] = 0.0; coords[1] = 0.0; coords[2] = 0.0; break; + case 3: coords[0] = 1.0; coords[1] = 0.0; coords[2] = 0.0; break; + case 6: coords[0] = 0.0; coords[1] = 0.5; coords[2] = 0.5; break; + case 5: coords[0] = 0.0; coords[1] = 0.0; coords[2] = 0.5; break; + case 4: coords[0] = 0.0; coords[1] = 0.5; coords[2] = 0.0; break; + case 7: coords[0] = 0.5; coords[1] = 0.5; coords[2] = 0.0; break; + case 9: coords[0] = 0.5; coords[1] = 0.0; coords[2] = 0.5; break; + case 8: coords[0] = 0.5; coords[1] = 0.0; coords[2] = 0.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = gc[1]*(2.0*gc[1] - 1.0); + funValue[2] = gc[2]*(2.0*gc[2] - 1.0); + funValue[1] = (1.0 - gc[0] - gc[1] - gc[2])*(1.0 - 2.0*gc[0] - 2.0*gc[1] - 2.0*gc[2]); + funValue[3] = gc[0]*(2.0*gc[0] - 1.0); + funValue[6] = 4.0*gc[1]*gc[2]; + funValue[5] = 4.0*gc[2]*(1.0 - gc[0] - gc[1] - gc[2]); + funValue[4] = 4.0*gc[1]*(1.0 - gc[0] - gc[1] - gc[2]); + funValue[7] = 4.0*gc[0]*gc[1]; + funValue[9] = 4.0*gc[0]*gc[2]; + funValue[8] = 4.0*gc[0]*(1.0 - gc[0] - gc[1] - gc[2]); + SHAPE_FUN_MACRO_END +} + +/*! + * Init Pyramid Reference coordinates ans Shape function. + * Case A. + */ +void GaussInfo::Pyra5aInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = 1.0; coords[1] = 0.0; coords[2] = 0.0; break; + case 1: coords[0] = 0.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 2: coords[0] = -1.0; coords[1] = 0.0; coords[2] = 0.0; break; + case 3: coords[0] = 0.0; coords[1] = -1.0; coords[2] = 0.0; break; + case 4: coords[0] = 0.0; coords[1] = 0.0; coords[2] = 1.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = 0.25*(-gc[0] + gc[1] - 1.0)*(-gc[0] - gc[1] - 1.0)*(1.0 - gc[2]); + funValue[1] = 0.25*(-gc[0] - gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]); + funValue[2] = 0.25*(+gc[0] + gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]); + funValue[3] = 0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0)*(1.0 - gc[2]); + funValue[4] = gc[2]; + SHAPE_FUN_MACRO_END +} +/*! + * Init Pyramid Reference coordinates ans Shape function. + * Case B. + */ +void GaussInfo::Pyra5bInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = 1.0; coords[1] = 0.0; coords[2] = 0.0; break; + case 3: coords[0] = 0.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 2: coords[0] = -1.0; coords[1] = 0.0; coords[2] = 0.0; break; + case 1: coords[0] = 0.0; coords[1] = -1.0; coords[2] = 0.0; break; + case 4: coords[0] = 0.0; coords[1] = 0.0; coords[2] = 1.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = 0.25*(-gc[0] + gc[1] - 1.0)*(-gc[0] - gc[1] - 1.0)*(1.0 - gc[2]); + funValue[3] = 0.25*(-gc[0] - gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]); + funValue[2] = 0.25*(+gc[0] + gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]); + funValue[1] = 0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0)*(1.0 - gc[2]); + funValue[4] = gc[2]; + SHAPE_FUN_MACRO_END +} + +/*! + * Init Quadratic Pyramid Reference coordinates ans Shape function. + * Case A. + */ +void GaussInfo::Pyra13aInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = 1.0; coords[1] = 0.0; coords[2] = 0.0; break; + case 1: coords[0] = 0.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 2: coords[0] = -1.0; coords[1] = 0.0; coords[2] = 0.0; break; + case 3: coords[0] = 0.0; coords[1] = -1.0; coords[2] = 0.0; break; + case 4: coords[0] = 0.0; coords[1] = 0.0; coords[2] = 1.0; break; + + case 5: coords[0] = 0.5; coords[1] = 0.5; coords[2] = 0.0; break; + case 6: coords[0] = -0.5; coords[1] = 0.5; coords[2] = 0.0; break; + case 7: coords[0] = -0.5; coords[1] = -0.5; coords[2] = 0.0; break; + case 8: coords[0] = 0.5; coords[1] = -0.5; coords[2] = 0.0; break; + case 9: coords[0] = 0.5; coords[1] = 0.0; coords[2] = 0.5; break; + case 10: coords[0] = 0.0; coords[1] = 0.5; coords[2] = 0.5; break; + case 11: coords[0] = -0.5; coords[1] = 0.0; coords[2] = 0.5; break; + case 12: coords[0] = 0.0; coords[1] = -0.5; coords[2] = 0.5; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)* + (gc[0] - 0.5)/(1.0 - gc[2]); + funValue[1] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] - gc[1] + gc[2] - 1.0)* + (gc[1] - 0.5)/(1.0 - gc[2]); + funValue[2] = 0.5*(+gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] + gc[1] + gc[2] - 1.0)* + (-gc[0] - 0.5)/(1.0 - gc[2]); + funValue[3] = 0.5*(+gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)* + (-gc[1] - 0.5)/(1.0 - gc[2]); + + funValue[4] = 2.0*gc[2]*(gc[2] - 0.5); + + funValue[5] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)* + (gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]); + funValue[6] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)* + (gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]); + funValue[7] = 0.5*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)* + (-gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]); + funValue[8] = 0.5*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)* + (-gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]); + + funValue[9] = 0.5*gc[2]*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)/ + (1.0 - gc[2]); + funValue[10] = 0.5*gc[2]*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)/ + (1.0 - gc[2]); + funValue[11] = 0.5*gc[2]*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)/ + (1.0 - gc[2]); + funValue[12] = 0.5*gc[2]*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)/ + (1.0 - gc[2]); + SHAPE_FUN_MACRO_END +} + +/*! + * Init Quadratic Pyramid Reference coordinates ans Shape function. + * Case B. + */ +void GaussInfo::Pyra13bInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = 1.0; coords[1] = 0.0; coords[2] = 0.0; break; + case 3: coords[0] = 0.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 2: coords[0] = -1.0; coords[1] = 0.0; coords[2] = 0.0; break; + case 1: coords[0] = 0.0; coords[1] = -1.0; coords[2] = 0.0; break; + case 4: coords[0] = 0.0; coords[1] = 0.0; coords[2] = 1.0; break; + case 8: coords[0] = 0.5; coords[1] = 0.5; coords[2] = 0.0; break; + case 7: coords[0] = -0.5; coords[1] = 0.5; coords[2] = 0.0; break; + case 6: coords[0] = -0.5; coords[1] = -0.5; coords[2] = 0.0; break; + case 5: coords[0] = 0.5; coords[1] = -0.5; coords[2] = 0.0; break; + case 9: coords[0] = 0.5; coords[1] = 0.0; coords[2] = 0.5; break; + case 12: coords[0] = 0.0; coords[1] = 0.5; coords[2] = 0.5; break; + case 11: coords[0] = -0.5; coords[1] = 0.0; coords[2] = 0.5; break; + case 10: coords[0] = 0.0; coords[1] = -0.5; coords[2] = 0.5; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)* + (gc[0] - 0.5)/(1.0 - gc[2]); + funValue[3] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] - gc[1] + gc[2] - 1.0)* + (gc[1] - 0.5)/(1.0 - gc[2]); + funValue[2] = 0.5*(+gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] + gc[1] + gc[2] - 1.0)* + (-gc[0] - 0.5)/(1.0 - gc[2]); + funValue[1] = 0.5*(+gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)* + (-gc[1] - 0.5)/(1.0 - gc[2]); + + funValue[4] = 2.0*gc[2]*(gc[2] - 0.5); + + funValue[8] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)* + (gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]); + funValue[7] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)* + (gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]); + funValue[6] = 0.5*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)* + (-gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]); + funValue[5] = 0.5*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)* + (-gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]); + + funValue[9] = 0.5*gc[2]*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)/ + (1.0 - gc[2]); + funValue[12] = 0.5*gc[2]*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)/ + (1.0 - gc[2]); + funValue[11] = 0.5*gc[2]*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)/ + (1.0 - gc[2]); + funValue[10] = 0.5*gc[2]*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)/ + (1.0 - gc[2]); + SHAPE_FUN_MACRO_END +} + + +/*! + * Init Pentahedron Reference coordinates and Shape function. + * Case A. + */ +void GaussInfo::Penta6aInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = -1.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 1: coords[0] = -1.0; coords[1] = -0.0; coords[2] = 1.0; break; + case 2: coords[0] = -1.0; coords[1] = 0.0; coords[2] = 0.0; break; + case 3: coords[0] = 1.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 4: coords[0] = 1.0; coords[1] = 0.0; coords[2] = 1.0; break; + case 5: coords[0] = 1.0; coords[1] = 0.0; coords[2] = 0.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = 0.5*gc[1]*(1.0 - gc[0]); + funValue[1] = 0.5*gc[2]*(1.0 - gc[0]); + funValue[2] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]); + + funValue[3] = 0.5*gc[1]*(gc[0] + 1.0); + funValue[4] = 0.5*gc[2]*(gc[0] + 1.0); + funValue[5] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]); + SHAPE_FUN_MACRO_END +} + +/*! + * Init Pentahedron Reference coordinates and Shape function. + * Case B. + */ +void GaussInfo::Penta6bInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = -1.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 2: coords[0] = -1.0; coords[1] = -0.0; coords[2] = 1.0; break; + case 1: coords[0] = -1.0; coords[1] = 0.0; coords[2] = 0.0; break; + case 3: coords[0] = 1.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 5: coords[0] = 1.0; coords[1] = 0.0; coords[2] = 1.0; break; + case 4: coords[0] = 1.0; coords[1] = 0.0; coords[2] = 0.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = 0.5*gc[1]*(1.0 - gc[0]); + funValue[2] = 0.5*gc[2]*(1.0 - gc[0]); + funValue[1] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]); + funValue[3] = 0.5*gc[1]*(gc[0] + 1.0); + funValue[5] = 0.5*gc[2]*(gc[0] + 1.0); + funValue[4] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]); + SHAPE_FUN_MACRO_END +} +/*! + * Init Pentahedron Reference coordinates and Shape function. + * Case A. + */ +void GaussInfo::Penta15aInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = -1.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 1: coords[0] = -1.0; coords[1] = -0.0; coords[2] = 1.0; break; + case 2: coords[0] = -1.0; coords[1] = 0.0; coords[2] = 0.0; break; + case 3: coords[0] = 1.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 4: coords[0] = 1.0; coords[1] = 0.0; coords[2] = 1.0; break; + case 5: coords[0] = 1.0; coords[1] = 0.0; coords[2] = 0.0; break; + + case 6: coords[0] = -1.0; coords[1] = 0.5; coords[2] = 0.5; break; + case 7: coords[0] = -1.0; coords[1] = 0.0; coords[2] = 0.5; break; + case 8: coords[0] = -1.0; coords[1] = 0.5; coords[2] = 0.0; break; + case 9: coords[0] = 0.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 10: coords[0] = 0.0; coords[1] = 0.0; coords[2] = 1.0; break; + case 11: coords[0] = 0.0; coords[1] = 0.0; coords[2] = 0.0; break; + case 12: coords[0] = 1.0; coords[1] = 0.5; coords[2] = 0.5; break; + case 13: coords[0] = 1.0; coords[1] = 0.0; coords[2] = 0.5; break; + case 14: coords[0] = 1.0; coords[1] = 0.5; coords[2] = 0.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]); + funValue[1] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]); + funValue[2] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]); + + funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]); + funValue[4] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]); + funValue[5] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]); + + funValue[6] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]); + funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]); + funValue[8] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]); + + funValue[9] = gc[1]*(1.0 - gc[0]*gc[0]); + funValue[10] = gc[2]*(1.0 - gc[0]*gc[0]); + funValue[11] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]); + + funValue[12] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]); + funValue[13] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]); + funValue[14] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]); + SHAPE_FUN_MACRO_END +} + +/*! + * Init Qaudratic Pentahedron Reference coordinates and Shape function. + * Case B. + */ +void GaussInfo::Penta15bInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = -1.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 2: coords[0] = -1.0; coords[1] = -0.0; coords[2] = 1.0; break; + case 1: coords[0] = -1.0; coords[1] = 0.0; coords[2] = 0.0; break; + case 3: coords[0] = 1.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 5: coords[0] = 1.0; coords[1] = 0.0; coords[2] = 1.0; break; + case 4: coords[0] = 1.0; coords[1] = 0.0; coords[2] = 0.0; break; + + case 8: coords[0] = -1.0; coords[1] = 0.5; coords[2] = 0.5; break; + case 7: coords[0] = -1.0; coords[1] = 0.0; coords[2] = 0.5; break; + case 6: coords[0] = -1.0; coords[1] = 0.5; coords[2] = 0.0; break; + case 12: coords[0] = 0.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 14: coords[0] = 0.0; coords[1] = 0.0; coords[2] = 1.0; break; + case 13: coords[0] = 0.0; coords[1] = 0.0; coords[2] = 0.0; break; + case 11: coords[0] = 1.0; coords[1] = 0.5; coords[2] = 0.5; break; + case 10: coords[0] = 1.0; coords[1] = 0.0; coords[2] = 0.5; break; + case 9: coords[0] = 1.0; coords[1] = 0.5; coords[2] = 0.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]); + funValue[2] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]); + funValue[1] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]); + + funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]); + funValue[5] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]); + funValue[4] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]); + + funValue[8] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]); + funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]); + funValue[6] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]); + + funValue[12] = gc[1]*(1.0 - gc[0]*gc[0]); + funValue[14] = gc[2]*(1.0 - gc[0]*gc[0]); + funValue[13] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]); + + funValue[11] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]); + funValue[10] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]); + funValue[9] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]); + SHAPE_FUN_MACRO_END +} + +/*! + * Init Hehahedron Reference coordinates and Shape function. + * Case A. + */ +void GaussInfo::Hexa8aInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = -1.0; coords[1] = -1.0; coords[2] = -1.0; break; + case 1: coords[0] = 1.0; coords[1] = -1.0; coords[2] = -1.0; break; + case 2: coords[0] = 1.0; coords[1] = 1.0; coords[2] = -1.0; break; + case 3: coords[0] = -1.0; coords[1] = 1.0; coords[2] = -1.0; break; + case 4: coords[0] = -1.0; coords[1] = -1.0; coords[2] = 1.0; break; + case 5: coords[0] = 1.0; coords[1] = -1.0; coords[2] = 1.0; break; + case 6: coords[0] = 1.0; coords[1] = 1.0; coords[2] = 1.0; break; + case 7: coords[0] = -1.0; coords[1] = 1.0; coords[2] = 1.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2]); + funValue[1] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2]); + funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2]); + funValue[3] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2]); + + funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2]); + funValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2]); + funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2]); + funValue[7] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2]); + SHAPE_FUN_MACRO_END +} + +/*! + * Init Hehahedron Reference coordinates and Shape function. + * Case B. + */ +void GaussInfo::Hexa8bInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = -1.0; coords[1] = -1.0; coords[2] = -1.0; break; + case 3: coords[0] = 1.0; coords[1] = -1.0; coords[2] = -1.0; break; + case 2: coords[0] = 1.0; coords[1] = 1.0; coords[2] = -1.0; break; + case 1: coords[0] = -1.0; coords[1] = 1.0; coords[2] = -1.0; break; + case 4: coords[0] = -1.0; coords[1] = -1.0; coords[2] = 1.0; break; + case 7: coords[0] = 1.0; coords[1] = -1.0; coords[2] = 1.0; break; + case 6: coords[0] = 1.0; coords[1] = 1.0; coords[2] = 1.0; break; + case 5: coords[0] = -1.0; coords[1] = 1.0; coords[2] = 1.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2]); + funValue[3] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2]); + funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2]); + funValue[1] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2]); + + funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2]); + funValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2]); + funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2]); + funValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2]); + SHAPE_FUN_MACRO_END +} + +/*! + * Init Qaudratic Hehahedron Reference coordinates and Shape function. + * Case A. + */ +void GaussInfo::Hexa20aInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = -1.0; coords[1] = -1.0; coords[2] = -1.0; break; + case 1: coords[0] = 1.0; coords[1] = -1.0; coords[2] = -1.0; break; + case 2: coords[0] = 1.0; coords[1] = 1.0; coords[2] = -1.0; break; + case 3: coords[0] = -1.0; coords[1] = 1.0; coords[2] = -1.0; break; + case 4: coords[0] = -1.0; coords[1] = -1.0; coords[2] = 1.0; break; + case 5: coords[0] = 1.0; coords[1] = -1.0; coords[2] = 1.0; break; + case 6: coords[0] = 1.0; coords[1] = 1.0; coords[2] = 1.0; break; + case 7: coords[0] = -1.0; coords[1] = 1.0; coords[2] = 1.0; break; + + case 8: coords[0] = 0.0; coords[1] = -1.0; coords[2] = -1.0; break; + case 9: coords[0] = 1.0; coords[1] = 0.0; coords[2] = -1.0; break; + case 10: coords[0] = 0.0; coords[1] = 1.0; coords[2] = -1.0; break; + case 11: coords[0] = -1.0; coords[1] = 0.0; coords[2] = -1.0; break; + case 12: coords[0] = -1.0; coords[1] = -1.0; coords[2] = 0.0; break; + case 13: coords[0] = 1.0; coords[1] = -1.0; coords[2] = 0.0; break; + case 14: coords[0] = 1.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 15: coords[0] = -1.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 16: coords[0] = 0.0; coords[1] = -1.0; coords[2] = 1.0; break; + case 17: coords[0] = 1.0; coords[1] = 0.0; coords[2] = 1.0; break; + case 18: coords[0] = 0.0; coords[1] = 1.0; coords[2] = 1.0; break; + case 19: coords[0] = -1.0; coords[1] = 0.0; coords[2] = 1.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2])* + (-2.0 - gc[0] - gc[1] - gc[2]); + funValue[1] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2])* + (-2.0 + gc[0] - gc[1] - gc[2]); + funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2])* + (-2.0 + gc[0] + gc[1] - gc[2]); + funValue[3] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2])* + (-2.0 - gc[0] + gc[1] - gc[2]); + funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2])* + (-2.0 - gc[0] - gc[1] + gc[2]); + funValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2])* + (-2.0 + gc[0] - gc[1] + gc[2]); + funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2])* + (-2.0 + gc[0] + gc[1] + gc[2]); + funValue[7] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2])* + (-2.0 - gc[0] + gc[1] + gc[2]); + + funValue[8] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]); + funValue[9] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]); + funValue[10] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]); + funValue[11] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]); + funValue[12] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]); + funValue[13] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]); + funValue[14] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]); + funValue[15] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]); + funValue[16] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]); + funValue[17] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]); + funValue[18] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]); + funValue[19] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]); + SHAPE_FUN_MACRO_END +} + +/*! + * Init Qaudratic Hehahedron Reference coordinates and Shape function. + * Case B. + */ +void GaussInfo::Hexa20bInit() { + LOCAL_COORD_MACRO_BEGIN + case 0: coords[0] = -1.0; coords[1] = -1.0; coords[2] = -1.0; break; + case 3: coords[0] = 1.0; coords[1] = -1.0; coords[2] = -1.0; break; + case 2: coords[0] = 1.0; coords[1] = 1.0; coords[2] = -1.0; break; + case 1: coords[0] = -1.0; coords[1] = 1.0; coords[2] = -1.0; break; + case 4: coords[0] = -1.0; coords[1] = -1.0; coords[2] = 1.0; break; + case 7: coords[0] = 1.0; coords[1] = -1.0; coords[2] = 1.0; break; + case 6: coords[0] = 1.0; coords[1] = 1.0; coords[2] = 1.0; break; + case 5: coords[0] = -1.0; coords[1] = 1.0; coords[2] = 1.0; break; + + case 11: coords[0] = 0.0; coords[1] = -1.0; coords[2] = -1.0; break; + case 10: coords[0] = 1.0; coords[1] = 0.0; coords[2] = -1.0; break; + case 9: coords[0] = 0.0; coords[1] = 1.0; coords[2] = -1.0; break; + case 8: coords[0] = -1.0; coords[1] = 0.0; coords[2] = -1.0; break; + case 16: coords[0] = -1.0; coords[1] = -1.0; coords[2] = 0.0; break; + case 19: coords[0] = 1.0; coords[1] = -1.0; coords[2] = 0.0; break; + case 18: coords[0] = 1.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 17: coords[0] = -1.0; coords[1] = 1.0; coords[2] = 0.0; break; + case 15: coords[0] = 0.0; coords[1] = -1.0; coords[2] = 1.0; break; + case 14: coords[0] = 1.0; coords[1] = 0.0; coords[2] = 1.0; break; + case 13: coords[0] = 0.0; coords[1] = 1.0; coords[2] = 1.0; break; + case 12: coords[0] = -1.0; coords[1] = 0.0; coords[2] = 1.0; break; + LOCAL_COORD_MACRO_END + + SHAPE_FUN_MACRO_BEGIN + + funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2])* + (-2.0 - gc[0] - gc[1] - gc[2]); + funValue[3] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2])* + (-2.0 + gc[0] - gc[1] - gc[2]); + funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2])* + (-2.0 + gc[0] + gc[1] - gc[2]); + funValue[1] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2])* + (-2.0 - gc[0] + gc[1] - gc[2]); + funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2])* + (-2.0 - gc[0] - gc[1] + gc[2]); + funValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2])* + (-2.0 + gc[0] - gc[1] + gc[2]); + funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2])* + (-2.0 + gc[0] + gc[1] + gc[2]); + funValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2])* + (-2.0 - gc[0] + gc[1] + gc[2]); + + funValue[11] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]); + funValue[10] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]); + funValue[9] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]); + funValue[8] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]); + funValue[16] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]); + funValue[19] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]); + funValue[18] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]); + funValue[17] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]); + funValue[15] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]); + funValue[14] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]); + funValue[13] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]); + funValue[12] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]); + SHAPE_FUN_MACRO_END +} + + + +//////////////////////////////////////////////////////////////////////////////////////////////// +// GAUSS COORD CLASS // +//////////////////////////////////////////////////////////////////////////////////////////////// +/*! + * Constructor + */ +GaussCoords::GaussCoords() { + +} + +/*! + * Destructor + */ +GaussCoords::~GaussCoords() { + GaussInfoVector::iterator it = myGaussInfo.begin(); + for( ; it != myGaussInfo.end(); it++ ) { + if((*it) != NULL) + delete (*it); + } +} + +/*! + * Add Gauss localization info + */ +void GaussCoords::addGaussInfo( NormalizedCellType theGeometry, + int coordDim, + const double* theGaussCoord, + int theNbGauss, + const double* theReferenceCoord, + int theNbRef) throw (INTERP_KERNEL::Exception) { + + GaussInfoVector::iterator it = myGaussInfo.begin(); + for( ; it != myGaussInfo.end(); it++ ) { + if( (*it)->GetCellType() == theGeometry ) { + break; + } + } + + DataVector aGaussCoord; + for(int i = 0 ; i < theNbGauss*coordDim; i++ ) + aGaussCoord.push_back(theGaussCoord[i]); + + DataVector aReferenceCoord; + for(int i = 0 ; i < theNbRef*coordDim; i++ ) + aReferenceCoord.push_back(theReferenceCoord[i]); + + + GaussInfo* info = new GaussInfo( theGeometry, aGaussCoord, theNbGauss, aReferenceCoord, theNbRef); + info->InitLocalInfo(); + + //If info with cell type doesn't exist add it + if( it == myGaussInfo.end() ) { + myGaussInfo.push_back(info); + + // If information exists, update it + } else { + int index = std::distance(myGaussInfo.begin(),it); + delete (*it); + myGaussInfo[index] = info; + } +} + + +/*! + * Calculate gauss points coordinates + */ +double* GaussCoords::CalculateCoords( NormalizedCellType theGeometry, + const double* theNodeCoords, + const int theSpaceDim, + const int* theIndex) + throw (INTERP_KERNEL::Exception) { + + GaussInfoVector::const_iterator it = myGaussInfo.begin(); + GaussInfoVector::const_iterator it_end = myGaussInfo.end(); + + //Try to find gauss localization info + for( ; it != it_end ; it++ ) { + if( (*it)->GetCellType() == theGeometry ) { + break; + } + } + if(it == it_end) { + throw INTERP_KERNEL::Exception("Can't find gauss localization information !!!"); + } + const GaussInfo* info = (*it); + int aConn = info->GetNbRef(); + + int nbCoords = theSpaceDim * info->GetNbGauss(); + double* aCoords = new double [nbCoords]; + for( int i = 0; i < nbCoords; i++ ) + aCoords[i] = 0.0; + +#ifdef MYDEBUG + std::cout<<"#######################################################"<GetCellType()<GetNbGauss() ; gaussId++ ) { +#ifdef MYDEBUG + std::cout<<"Gauss ID = "<GetFunctionValues(gaussId); + for ( int connId = 0; connId < aConn ; connId++ ) { + const double* nodeCoord = theNodeCoords + (theIndex[connId]*theSpaceDim); + +#ifdef MYDEBUG + std::cout<<"Function Value["< + +namespace INTERP_KERNEL { + + typedef std::vector DataVector; + typedef std::vector IndexVector; + + //Class to store Gauss Points information + class GaussInfo { + + public: + GaussInfo( NormalizedCellType theGeometry, + const DataVector& theGaussCoord, + int theNbGauss, + const DataVector& theReferenceCoord, + int theNbRef + ); + ~GaussInfo(); + + NormalizedCellType GetCellType() const; + + int GetGaussCoordDim() const; + int GetReferenceCoordDim() const; + + int GetNbGauss() const; + int GetNbRef() const; + + const double* GetFunctionValues( const int theGaussId ) const; + + void InitLocalInfo() throw (INTERP_KERNEL::Exception); + + protected: + + bool isSatisfy(); + + //1D + void Seg2Init(); + void Seg3Init(); + + //2D + void Tria3aInit(); + void Tria3bInit(); + void Tria6aInit(); + void Tria6bInit(); + + void Quad4aInit(); + void Quad4bInit(); + void Quad8aInit(); + void Quad8bInit(); + + //3D + void Tetra4aInit(); + void Tetra4bInit(); + void Tetra10aInit(); + void Tetra10bInit(); + + void Pyra5aInit(); + void Pyra5bInit(); + void Pyra13aInit(); + void Pyra13bInit(); + + void Penta6aInit(); + void Penta6bInit(); + void Penta15aInit(); + void Penta15bInit(); + + void Hexa8aInit(); + void Hexa8bInit(); + void Hexa20aInit(); + void Hexa20bInit(); + + + private: + //INFORMATION from MEDMEM + NormalizedCellType myGeometry; //Cell type + + int myNbGauss; //Nb of the gauss points for element + DataVector myGaussCoord; //Gauss coordinates + + int myNbRef; //Nb of the nodes for element: + //NORM_SEG2 - 2 + //NORM_SEG3 - 3 + //NORM_TRI3 - 3 + //............. + + DataVector myReferenceCoord; //Reference coordinates + + //LOCAL INFORMATION + DataVector myLocalReferenceCoord; //Vector to store reference coordinates + int myLocalRefDim; //Dimension of the local reference coordinates: + // (x) - 1D case + // (x, y) - 2D case + // (x, y, z) - 3D case + int myLocalNbRef; //Nb of the local reference coordinates + + DataVector myFunctionValue; //Shape Function values + }; + + + //Class for calculation of the coordinates of the gauss points + class GaussCoords { + public: + + GaussCoords(); + ~GaussCoords(); + + void addGaussInfo( NormalizedCellType theGeometry, + int coordDim, + const double* theGaussCoord, + int theNbGauss, + const double* theReferenceCoord, + int theNbRef) throw (INTERP_KERNEL::Exception); + + double* CalculateCoords( NormalizedCellType theGeometry, + const double* theNodeCoords, + const int theDimSpace, + const int* theIndex + ) throw(INTERP_KERNEL::Exception); + private: + typedef std::vector GaussInfoVector; + GaussInfoVector myGaussInfo; + }; +} +#endif //INTERPKERNELGAUSS diff --git a/src/INTERP_KERNEL/GaussPoints/Makefile.am b/src/INTERP_KERNEL/GaussPoints/Makefile.am new file mode 100644 index 000000000..0a95f62c3 --- /dev/null +++ b/src/INTERP_KERNEL/GaussPoints/Makefile.am @@ -0,0 +1,37 @@ +# Copyright (C) 2007-2010 CEA/DEN, EDF R&D +# +# 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. +# +# 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 +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# + +# File : Makefile.am +# Author : Vincent BERGEAUD (CEA/DEN/DANS/DM2S/SFME/LGLS) +# Module : MED +# +include $(top_srcdir)/adm_local/unix/make_common_starter.am + +noinst_LTLIBRARIES = libinterpkernelgauss.la + +salomeinclude_HEADERS = \ +GaussCoords.hxx + +# Libraries targets +dist_libinterpkernelgauss_la_SOURCES = \ +GaussCoords.cxx + +libinterpkernelgauss_la_CPPFLAGS=-I$(srcdir)/../Bases -I$(srcdir)/.. + +AM_CPPFLAGS += $(libinterpkernelgauss_la_CPPFLAGS) \ No newline at end of file diff --git a/src/INTERP_KERNEL/Makefile.am b/src/INTERP_KERNEL/Makefile.am index 53a6650f0..4268d2c61 100644 --- a/src/INTERP_KERNEL/Makefile.am +++ b/src/INTERP_KERNEL/Makefile.am @@ -24,9 +24,9 @@ # include $(top_srcdir)/adm_local/unix/make_common_starter.am -SUBDIRS = Bases Geometric2D ExprEval . +SUBDIRS = Bases Geometric2D ExprEval GaussPoints . -DIST_SUBDIRS = Bases Geometric2D ExprEval +DIST_SUBDIRS = Bases Geometric2D ExprEval GaussPoints lib_LTLIBRARIES = libinterpkernel.la @@ -198,7 +198,8 @@ libinterpkernel_la_CPPFLAGS=-I$(srcdir)/Geometric2D -I$(srcdir)/Bases libinterpkernel_la_LDFLAGS= # the geom2D library is included in the interpkernel one -libinterpkernel_la_LIBADD= ./Geometric2D/libInterpGeometric2DAlg.la Bases/libinterpkernelbases.la ExprEval/libinterpkernelexpreval.la +libinterpkernel_la_LIBADD= ./Geometric2D/libInterpGeometric2DAlg.la Bases/libinterpkernelbases.la ExprEval/libinterpkernelexpreval.la \ + GaussPoints/libinterpkernelgauss.la AM_CPPFLAGS += $(libinterpkernel_la_CPPFLAGS) LDADD= $(libinterpkernel_la_LDFLAGS) diff --git a/src/INTERP_KERNELTest/Makefile.am b/src/INTERP_KERNELTest/Makefile.am index 12a36ef3e..61fcfca0d 100644 --- a/src/INTERP_KERNELTest/Makefile.am +++ b/src/INTERP_KERNELTest/Makefile.am @@ -67,6 +67,7 @@ libInterpKernelTest_la_CPPFLAGS = \ -I$(srcdir)/../INTERP_KERNEL/Geometric2D \ -I$(srcdir)/../INTERP_KERNEL/Bases \ -I$(srcdir)/../INTERP_KERNEL/ExprEval \ + -I$(srcdir)/../INTERP_KERNEL/GaussPoints \ -DOPTIMIZE -DLOG_LEVEL=0 libInterpKernelTest_la_LDFLAGS = \ diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx index 72df7f58c..eb829e58e 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2010 CEA/DEN, EDF R&D + // Copyright (C) 2007-2010 CEA/DEN, EDF R&D // // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public diff --git a/src/MEDCoupling_Swig/MEDCoupling.i b/src/MEDCoupling_Swig/MEDCoupling.i index 94c2613d4..87b1562c9 100644 --- a/src/MEDCoupling_Swig/MEDCoupling.i +++ b/src/MEDCoupling_Swig/MEDCoupling.i @@ -66,8 +66,8 @@ using namespace INTERP_KERNEL; %newobject ParaMEDMEM::MEDCouplingFieldDouble::New; %newobject ParaMEDMEM::MEDCouplingFieldDouble::getArray; %newobject ParaMEDMEM::MEDCouplingFieldDouble::getEndArray; -%newobject ParaMEDMEM::MEDCouplingFieldDouble::mergeFields; -%newobject ParaMEDMEM::MEDCouplingFieldDouble::meldFields; +%newobject ParaMEDMEM::MEDCouplingFieldDouble::MergeFields; +%newobject ParaMEDMEM::MEDCouplingFieldDouble::MeldFields; %newobject ParaMEDMEM::MEDCouplingFieldDouble::doublyContractedProduct; %newobject ParaMEDMEM::MEDCouplingFieldDouble::determinant; %newobject ParaMEDMEM::MEDCouplingFieldDouble::eigenValues; @@ -78,13 +78,17 @@ using namespace INTERP_KERNEL; %newobject ParaMEDMEM::MEDCouplingFieldDouble::magnitude; %newobject ParaMEDMEM::MEDCouplingFieldDouble::maxPerTuple; %newobject ParaMEDMEM::MEDCouplingFieldDouble::keepSelectedComponents; -%newobject ParaMEDMEM::MEDCouplingFieldDouble::dotFields; +%newobject ParaMEDMEM::MEDCouplingFieldDouble::DotFields; %newobject ParaMEDMEM::MEDCouplingFieldDouble::dot; -%newobject ParaMEDMEM::MEDCouplingFieldDouble::crossProductFields; +%newobject ParaMEDMEM::MEDCouplingFieldDouble::CrossProductFields; %newobject ParaMEDMEM::MEDCouplingFieldDouble::crossProduct; -%newobject ParaMEDMEM::MEDCouplingFieldDouble::maxFields; +%newobject ParaMEDMEM::MEDCouplingFieldDouble::MaxFields; %newobject ParaMEDMEM::MEDCouplingFieldDouble::max; -%newobject ParaMEDMEM::MEDCouplingFieldDouble::minFields; +%newobject ParaMEDMEM::MEDCouplingFieldDouble::MinFields; +%newobject ParaMEDMEM::MEDCouplingFieldDouble::AddFields; +%newobject ParaMEDMEM::MEDCouplingFieldDouble::SubstractFields; +%newobject ParaMEDMEM::MEDCouplingFieldDouble::MultiplyFields; +%newobject ParaMEDMEM::MEDCouplingFieldDouble::DivideFields; %newobject ParaMEDMEM::MEDCouplingFieldDouble::min; %newobject ParaMEDMEM::MEDCouplingFieldDouble::getIdsInRange; %newobject ParaMEDMEM::MEDCouplingFieldDouble::buildSubPart; @@ -112,8 +116,10 @@ using namespace INTERP_KERNEL; %newobject ParaMEDMEM::DataArrayInt::invertArrayN2O2O2N; %newobject ParaMEDMEM::DataArrayInt::getIdsEqual; %newobject ParaMEDMEM::DataArrayInt::getIdsEqualList; -%newobject ParaMEDMEM::DataArrayInt::aggregate; -%newobject ParaMEDMEM::DataArrayInt::meld; +%newobject ParaMEDMEM::DataArrayInt::Aggregate; +%newobject ParaMEDMEM::DataArrayInt::Meld; +%newobject ParaMEDMEM::DataArrayInt::BuildUnion; +%newobject ParaMEDMEM::DataArrayInt::BuildIntersection; %newobject ParaMEDMEM::DataArrayInt::fromNoInterlace; %newobject ParaMEDMEM::DataArrayInt::toNoInterlace; %newobject ParaMEDMEM::DataArrayInt::buildComplement; @@ -121,18 +127,19 @@ using namespace INTERP_KERNEL; %newobject ParaMEDMEM::DataArrayInt::buildSubstraction; %newobject ParaMEDMEM::DataArrayInt::buildIntersection; %newobject ParaMEDMEM::DataArrayInt::deltaShiftIndex; +%newobject ParaMEDMEM::DataArrayInt::buildPermutationArr; %newobject ParaMEDMEM::DataArrayDouble::New; %newobject ParaMEDMEM::DataArrayDouble::convertToIntArr; %newobject ParaMEDMEM::DataArrayDouble::deepCpy; %newobject ParaMEDMEM::DataArrayDouble::performCpy; -%newobject ParaMEDMEM::DataArrayDouble::aggregate; -%newobject ParaMEDMEM::DataArrayDouble::meld; -%newobject ParaMEDMEM::DataArrayDouble::dot; -%newobject ParaMEDMEM::DataArrayDouble::crossProduct; -%newobject ParaMEDMEM::DataArrayDouble::add; -%newobject ParaMEDMEM::DataArrayDouble::substract; -%newobject ParaMEDMEM::DataArrayDouble::multiply; -%newobject ParaMEDMEM::DataArrayDouble::divide; +%newobject ParaMEDMEM::DataArrayDouble::Aggregate; +%newobject ParaMEDMEM::DataArrayDouble::Meld; +%newobject ParaMEDMEM::DataArrayDouble::Dot; +%newobject ParaMEDMEM::DataArrayDouble::CrossProduct; +%newobject ParaMEDMEM::DataArrayDouble::Add; +%newobject ParaMEDMEM::DataArrayDouble::Substract; +%newobject ParaMEDMEM::DataArrayDouble::Multiply; +%newobject ParaMEDMEM::DataArrayDouble::Divide; %newobject ParaMEDMEM::DataArrayDouble::substr; %newobject ParaMEDMEM::DataArrayDouble::changeNbOfComponents; %newobject ParaMEDMEM::DataArrayDouble::keepSelectedComponents; @@ -166,19 +173,20 @@ using namespace INTERP_KERNEL; %newobject ParaMEDMEM::MEDCouplingMesh::getMeasureField; %newobject ParaMEDMEM::MEDCouplingMesh::simplexize; %newobject ParaMEDMEM::MEDCouplingMesh::buildUnstructured; +%newobject ParaMEDMEM::MEDCouplingMesh::MergeMeshes; %newobject ParaMEDMEM::MEDCouplingPointSet::zipCoordsTraducer; %newobject ParaMEDMEM::MEDCouplingPointSet::buildBoundaryMesh; -%newobject ParaMEDMEM::MEDCouplingPointSet::mergeNodesArray; -%newobject ParaMEDMEM::MEDCouplingPointSet::buildInstanceFromMeshType; +%newobject ParaMEDMEM::MEDCouplingPointSet::MergeNodesArray; +%newobject ParaMEDMEM::MEDCouplingPointSet::BuildInstanceFromMeshType; %newobject ParaMEDMEM::MEDCouplingUMesh::New; %newobject ParaMEDMEM::MEDCouplingUMesh::getNodalConnectivity; %newobject ParaMEDMEM::MEDCouplingUMesh::getNodalConnectivityIndex; %newobject ParaMEDMEM::MEDCouplingUMesh::clone; %newobject ParaMEDMEM::MEDCouplingUMesh::zipConnectivityTraducer; %newobject ParaMEDMEM::MEDCouplingUMesh::buildDescendingConnectivity; -%newobject ParaMEDMEM::MEDCouplingUMesh::buildExtrudedMeshFromThis; -%newobject ParaMEDMEM::MEDCouplingUMesh::mergeUMeshes; -%newobject ParaMEDMEM::MEDCouplingUMesh::mergeUMeshesOnSameCoords; +%newobject ParaMEDMEM::MEDCouplingUMesh::buildExtrudedMesh; +%newobject ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes; +%newobject ParaMEDMEM::MEDCouplingUMesh::MergeUMeshesOnSameCoords; %newobject ParaMEDMEM::MEDCouplingUMesh::buildNewNumberingFromCommNodesFrmt; %newobject ParaMEDMEM::MEDCouplingUMesh::rearrange2ConsecutiveCellTypes; %newobject ParaMEDMEM::MEDCouplingUMesh::convertCellArrayPerGeoType; @@ -273,7 +281,7 @@ namespace ParaMEDMEM virtual MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const throw(INTERP_KERNEL::Exception) = 0; virtual bool areCompatibleForMerge(const MEDCouplingMesh *other) const throw(INTERP_KERNEL::Exception); virtual DataArrayInt *simplexize(int policy) throw(INTERP_KERNEL::Exception); - static MEDCouplingMesh *mergeMeshes(const MEDCouplingMesh *mesh1, const MEDCouplingMesh *mesh2) throw(INTERP_KERNEL::Exception); + static MEDCouplingMesh *MergeMeshes(const MEDCouplingMesh *mesh1, const MEDCouplingMesh *mesh2) throw(INTERP_KERNEL::Exception); %extend { std::string __str__() const @@ -409,8 +417,8 @@ namespace ParaMEDMEM void changeSpaceDimension(int newSpaceDim, double dftVal=0.) throw(INTERP_KERNEL::Exception); void tryToShareSameCoords(const MEDCouplingPointSet& other, double epsilon) throw(INTERP_KERNEL::Exception); virtual void tryToShareSameCoordsPermute(const MEDCouplingPointSet& other, double epsilon) throw(INTERP_KERNEL::Exception) = 0; - static DataArrayDouble *mergeNodesArray(const MEDCouplingPointSet *m1, const MEDCouplingPointSet *m2) throw(INTERP_KERNEL::Exception); - static MEDCouplingPointSet *buildInstanceFromMeshType(MEDCouplingMeshType type) throw(INTERP_KERNEL::Exception); + static DataArrayDouble *MergeNodesArray(const MEDCouplingPointSet *m1, const MEDCouplingPointSet *m2) throw(INTERP_KERNEL::Exception); + static MEDCouplingPointSet *BuildInstanceFromMeshType(MEDCouplingMeshType type) throw(INTERP_KERNEL::Exception); virtual MEDCouplingPointSet *buildBoundaryMesh(bool keepCoords) const throw(INTERP_KERNEL::Exception) = 0; virtual bool isEmptyMesh(const std::vector& tinyInfo) const throw(INTERP_KERNEL::Exception) = 0; //! size of returned tinyInfo must be always the same. @@ -563,24 +571,24 @@ namespace ParaMEDMEM return convertIntArrToPyList2(elems); } - static void rotate2DAlg(PyObject *center, double angle, int nbNodes, PyObject *coords) throw(INTERP_KERNEL::Exception) + static void Rotate2DAlg(PyObject *center, double angle, int nbNodes, PyObject *coords) throw(INTERP_KERNEL::Exception) { int sz; double *c=convertPyToNewDblArr2(center,&sz); double *coo=convertPyToNewDblArr2(coords,&sz); - ParaMEDMEM::MEDCouplingPointSet::rotate2DAlg(c,angle,nbNodes,coo); + ParaMEDMEM::MEDCouplingPointSet::Rotate2DAlg(c,angle,nbNodes,coo); for(int i=0;i meshes; - if(PyList_Check(ms)) - { - int sz=PyList_Size(ms); - meshes.resize(sz); - for(int i=0;i(arg); - } - } - else - { - PyErr_SetString(PyExc_TypeError,"mergeUMeshesOnSameCoords : not a list as first parameter"); - PyErr_Print(); - return 0; - } - MEDCouplingUMesh *ret=MEDCouplingUMesh::mergeUMeshesOnSameCoords(meshes); + std::vector meshes; + convertPyObjToVecUMeshesCst(ms,meshes); + MEDCouplingUMesh *ret=MEDCouplingUMesh::MergeUMeshesOnSameCoords(meshes); return convertMesh(ret, SWIG_POINTER_OWN | 0 ); } - static PyObject *fuseUMeshesOnSameCoords(PyObject *ms, int compType) throw(INTERP_KERNEL::Exception) + static PyObject *FuseUMeshesOnSameCoords(PyObject *ms, int compType) throw(INTERP_KERNEL::Exception) { int sz; - std::vector meshes; - convertPyObjToVecUMeshes(ms,meshes); + std::vector meshes; + convertPyObjToVecUMeshesCst(ms,meshes); std::vector corr; - MEDCouplingUMesh *um=MEDCouplingUMesh::fuseUMeshesOnSameCoords(meshes,compType,corr); + MEDCouplingUMesh *um=MEDCouplingUMesh::FuseUMeshesOnSameCoords(meshes,compType,corr); sz=corr.size(); PyObject *ret1=PyList_New(sz); for(int i=0;i tmp; convertPyObjToVecUMeshesCst(li,tmp); - return MEDCouplingUMesh::mergeUMeshes(tmp); + return MEDCouplingUMesh::MergeUMeshes(tmp); } PyObject *areCellsIncludedIn(const MEDCouplingUMesh *other, int compType) const throw(INTERP_KERNEL::Exception) @@ -881,7 +872,7 @@ namespace ParaMEDMEM } void convertToPolyTypes(const std::vector& cellIdsToConvert) throw(INTERP_KERNEL::Exception); void unPolyze() throw(INTERP_KERNEL::Exception); - MEDCouplingUMesh *buildExtrudedMeshFromThis(const MEDCouplingUMesh *mesh1D, int policy) throw(INTERP_KERNEL::Exception); + MEDCouplingUMesh *buildExtrudedMesh(const MEDCouplingUMesh *mesh1D, int policy) throw(INTERP_KERNEL::Exception); }; class MEDCouplingExtrudedMesh : public ParaMEDMEM::MEDCouplingMesh @@ -965,8 +956,16 @@ namespace ParaMEDMEM void setValues(PyObject *li, int nbOfTuples, int nbOfElsPerTuple) throw(INTERP_KERNEL::Exception) { - int sz; - double *tmp=convertPyToNewDblArr2(li,&sz); + double *tmp=new double[nbOfTuples*nbOfElsPerTuple]; + try + { + fillArrayWithPyListDbl(li,tmp,nbOfTuples*nbOfElsPerTuple,0.); + } + catch(INTERP_KERNEL::Exception& e) + { + delete [] tmp; + throw e; + } self->useArray(tmp,true,CPP_DEALLOC,nbOfTuples,nbOfElsPerTuple); } @@ -1138,18 +1137,18 @@ namespace ParaMEDMEM return convertDblArrToPyList(tmp,sz); } - static DataArrayDouble *aggregate(PyObject *li) throw(INTERP_KERNEL::Exception) + static DataArrayDouble *Aggregate(PyObject *li) throw(INTERP_KERNEL::Exception) { std::vector tmp; convertPyObjToVecDataArrayDblCst(li,tmp); - return DataArrayDouble::aggregate(tmp); + return DataArrayDouble::Aggregate(tmp); } - static DataArrayDouble *meld(PyObject *li) throw(INTERP_KERNEL::Exception) + static DataArrayDouble *Meld(PyObject *li) throw(INTERP_KERNEL::Exception) { std::vector tmp; convertPyObjToVecDataArrayDblCst(li,tmp); - return DataArrayDouble::meld(tmp); + return DataArrayDouble::Meld(tmp); } }; @@ -1160,10 +1159,24 @@ namespace ParaMEDMEM return self->repr(); } + PyObject *getDifferentValues(bool val) const throw(INTERP_KERNEL::Exception) + { + std::set ret=self->getDifferentValues(); + return convertIntArrToPyList3(ret); + } + void setValues(PyObject *li, int nbOfTuples, int nbOfElsPerTuple) throw(INTERP_KERNEL::Exception) { - int size; - int *tmp=convertPyToNewIntArr2(li,&size); + int *tmp=new int[nbOfTuples*nbOfElsPerTuple]; + try + { + fillArrayWithPyListInt(li,tmp,nbOfTuples*nbOfElsPerTuple,0.); + } + catch(INTERP_KERNEL::Exception& e) + { + delete [] tmp; + throw e; + } self->useArray(tmp,true,CPP_DEALLOC,nbOfTuples,nbOfElsPerTuple); } @@ -1181,12 +1194,12 @@ namespace ParaMEDMEM return convertIntArrToPyListOfTuple(vals,nbOfComp,nbOfTuples); } - static PyObject *makePartition(PyObject *gps, int newNb) throw(INTERP_KERNEL::Exception) + static PyObject *MakePartition(PyObject *gps, int newNb) throw(INTERP_KERNEL::Exception) { - std::vector groups; + std::vector groups; std::vector< std::vector > fidsOfGroups; - convertPyObjToVecDataArrayInt(gps,groups); - ParaMEDMEM::DataArrayInt *ret0=ParaMEDMEM::DataArrayInt::makePartition(groups,newNb,fidsOfGroups); + convertPyObjToVecDataArrayIntCst(gps,groups); + ParaMEDMEM::DataArrayInt *ret0=ParaMEDMEM::DataArrayInt::MakePartition(groups,newNb,fidsOfGroups); PyObject *ret = PyList_New(2); PyList_SetItem(ret,0,SWIG_NewPointerObj(SWIG_as_voidptr(ret0),SWIGTYPE_p_ParaMEDMEM__DataArrayInt, SWIG_POINTER_OWN | 0 )); int sz=fidsOfGroups.size(); @@ -1293,11 +1306,25 @@ namespace ParaMEDMEM return convertIntArrToPyList(tmp,sz); } - static DataArrayInt *meld(PyObject *li) throw(INTERP_KERNEL::Exception) + static DataArrayInt *Meld(PyObject *li) throw(INTERP_KERNEL::Exception) + { + std::vector tmp; + convertPyObjToVecDataArrayIntCst(li,tmp); + return DataArrayInt::Meld(tmp); + } + + static DataArrayInt *BuildUnion(PyObject *li) throw(INTERP_KERNEL::Exception) + { + std::vector tmp; + convertPyObjToVecDataArrayIntCst(li,tmp); + return DataArrayInt::BuildUnion(tmp); + } + + static DataArrayInt *BuildIntersection(PyObject *li) throw(INTERP_KERNEL::Exception) { std::vector tmp; convertPyObjToVecDataArrayIntCst(li,tmp); - return DataArrayInt::meld(tmp); + return DataArrayInt::BuildIntersection(tmp); } PyObject *getMaxValue() const throw(INTERP_KERNEL::Exception) @@ -1417,6 +1444,12 @@ namespace ParaMEDMEM void setNature(NatureOfField nat) throw(INTERP_KERNEL::Exception); void setTimeTolerance(double val) throw(INTERP_KERNEL::Exception); double getTimeTolerance() const throw(INTERP_KERNEL::Exception); + void setIteration(int it) throw(INTERP_KERNEL::Exception); + void setEndIteration(int it) throw(INTERP_KERNEL::Exception); + void setOrder(int order) throw(INTERP_KERNEL::Exception); + void setEndOrder(int order) throw(INTERP_KERNEL::Exception); + void setTimeValue(double val) throw(INTERP_KERNEL::Exception); + void setEndTimeValue(double val) throw(INTERP_KERNEL::Exception); void updateTime() throw(INTERP_KERNEL::Exception); void changeUnderlyingMesh(const MEDCouplingMesh *other, int levOfCheck, double prec) throw(INTERP_KERNEL::Exception); void substractInPlaceDM(const MEDCouplingFieldDouble *f, int levOfCheck, double prec) throw(INTERP_KERNEL::Exception); @@ -1453,15 +1486,19 @@ namespace ParaMEDMEM double normL2(int compId) const throw(INTERP_KERNEL::Exception); DataArrayInt *getIdsInRange(double vmin, double vmax) const throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *buildSubPart(const DataArrayInt *part) const throw(INTERP_KERNEL::Exception); - static MEDCouplingFieldDouble *mergeFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception); - static MEDCouplingFieldDouble *meldFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception); - static MEDCouplingFieldDouble *dotFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception); + static MEDCouplingFieldDouble *MergeFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception); + static MEDCouplingFieldDouble *MeldFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception); + static MEDCouplingFieldDouble *DotFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *dot(const MEDCouplingFieldDouble& other) const throw(INTERP_KERNEL::Exception); - static MEDCouplingFieldDouble *crossProductFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception); + static MEDCouplingFieldDouble *CrossProductFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *crossProduct(const MEDCouplingFieldDouble& other) const throw(INTERP_KERNEL::Exception); - static MEDCouplingFieldDouble *maxFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception); + static MEDCouplingFieldDouble *MaxFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *max(const MEDCouplingFieldDouble& other) const throw(INTERP_KERNEL::Exception); - static MEDCouplingFieldDouble *minFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception); + static MEDCouplingFieldDouble *MinFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception); + static MEDCouplingFieldDouble *AddFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception); + static MEDCouplingFieldDouble *SubstractFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception); + static MEDCouplingFieldDouble *MultiplyFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception); + static MEDCouplingFieldDouble *DivideFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *min(const MEDCouplingFieldDouble& other) const throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *operator+(const MEDCouplingFieldDouble& other) const throw(INTERP_KERNEL::Exception); const MEDCouplingFieldDouble &operator+=(const MEDCouplingFieldDouble& other) throw(INTERP_KERNEL::Exception); @@ -1694,11 +1731,11 @@ namespace ParaMEDMEM self->setSelectedComponents(f,tmp); } - static MEDCouplingFieldDouble *mergeFields(PyObject *li) throw(INTERP_KERNEL::Exception) + static MEDCouplingFieldDouble *MergeFields(PyObject *li) throw(INTERP_KERNEL::Exception) { std::vector tmp; convertPyObjToVecFieldDblCst(li,tmp); - return MEDCouplingFieldDouble::mergeFields(tmp); + return MEDCouplingFieldDouble::MergeFields(tmp); } } }; diff --git a/src/MEDCoupling_Swig/MEDCouplingTypemaps.i b/src/MEDCoupling_Swig/MEDCouplingTypemaps.i index 04254c96d..a42063aec 100644 --- a/src/MEDCoupling_Swig/MEDCouplingTypemaps.i +++ b/src/MEDCoupling_Swig/MEDCouplingTypemaps.i @@ -70,6 +70,16 @@ static PyObject *convertIntArrToPyList2(const std::vector& v) throw(INTERP_ #endif } +static PyObject *convertIntArrToPyList3(const std::set& v) throw(INTERP_KERNEL::Exception) +{ + int size=v.size(); + PyObject *ret=PyList_New(size); + std::set::const_iterator it=v.begin(); + for(int i=0;i& arr) throw(I } } +static void fillArrayWithPyListInt(PyObject *pyLi, int *arrToFill, int sizeOfArray, int dftVal) throw(INTERP_KERNEL::Exception) +{ + if(PyList_Check(pyLi)) + { + int size=PyList_Size(pyLi); + for(int i=0;i& v) throw(INTERP_KERNEL::Exception) +static void fillArrayWithPyListDbl(PyObject *pyLi, double *arrToFill, int sizeOfArray, double dftVal) throw(INTERP_KERNEL::Exception) { - if(PyList_Check(ms)) + if(PyList_Check(pyLi)) { - int size=PyList_Size(ms); - v.resize(size); + int size=PyList_Size(pyLi); for(int i=0;i(argp); - v[i]=arg; } + for(int i=size;i& v) throw(INTERP_KERNEL::Exception) -{ - if(PyList_Check(ms)) - { - int size=PyList_Size(ms); - v.resize(size); - for(int i=0;i(argp); - v[i]=arg; - } - } - else - { - const char msg[]="convertPyObjToVecDataArrayInt : not a list"; - PyErr_SetString(PyExc_TypeError,msg); - throw INTERP_KERNEL::Exception(msg); - } -} - void convertPyObjToVecDataArrayIntCst(PyObject *ms, std::vector& v) throw(INTERP_KERNEL::Exception) { if(PyList_Check(ms)) diff --git a/src/MEDLoader/MEDFileMesh.cxx b/src/MEDLoader/MEDFileMesh.cxx index 89f93a3d3..5499fcbd7 100644 --- a/src/MEDLoader/MEDFileMesh.cxx +++ b/src/MEDLoader/MEDFileMesh.cxx @@ -736,6 +736,22 @@ void MEDFileUMesh::addNodeGroup(const std::string& name, const std::vector& } +void MEDFileUMesh::copyFamGrpMapsFrom(const MEDFileUMesh& other) +{ + _groups=other._groups; + _families=other._families; +} + +void MEDFileUMesh::setFamilyInfo(const std::map& info) +{ + _families=info; +} + +void MEDFileUMesh::setGroupInfo(const std::map >&info) +{ + _groups=info; +} + void MEDFileUMesh::setFamilyNameAttachedOnId(int id, const std::string& newFamName) throw(INTERP_KERNEL::Exception) { std::string oldName=getFamilyNameGivenId(id); diff --git a/src/MEDLoader/MEDFileMesh.hxx b/src/MEDLoader/MEDFileMesh.hxx index 1e99185af..cd103ff92 100644 --- a/src/MEDLoader/MEDFileMesh.hxx +++ b/src/MEDLoader/MEDFileMesh.hxx @@ -91,9 +91,14 @@ namespace ParaMEDMEM MEDCouplingUMesh *getLevelM1Mesh(bool renum=true) const throw(INTERP_KERNEL::Exception); MEDCouplingUMesh *getLevelM2Mesh(bool renum=true) const throw(INTERP_KERNEL::Exception); MEDCouplingUMesh *getLevelM3Mesh(bool renum=true) const throw(INTERP_KERNEL::Exception); + const std::map& getFamilyInfo() const { return _families; } + const std::map >& getGroupInfo() const { return _groups; } bool existsFamily(int famId) const; bool existsFamily(const char *familyName) const; // + void copyFamGrpMapsFrom(const MEDFileUMesh& other); + void setFamilyInfo(const std::map& info); + void setGroupInfo(const std::map >&info); void setFamilyNameAttachedOnId(int id, const std::string& newFamName) throw(INTERP_KERNEL::Exception); void setCoords(DataArrayDouble *coords) throw(INTERP_KERNEL::Exception); void setGroupsAtLevel(int meshDimRelToMaxExt, const std::vector& grps, bool renum=true) throw(INTERP_KERNEL::Exception); diff --git a/src/MEDLoader/Swig/MEDLoader.i b/src/MEDLoader/Swig/MEDLoader.i index e6f960b7a..10b8168ad 100644 --- a/src/MEDLoader/Swig/MEDLoader.i +++ b/src/MEDLoader/Swig/MEDLoader.i @@ -63,11 +63,11 @@ using namespace ParaMEDMEM; %newobject ParaMEDMEM::MEDFileUMesh::getNodeGroupsArr; %newobject ParaMEDMEM::MEDFileUMesh::getNodeFamilyArr; %newobject ParaMEDMEM::MEDFileUMesh::getNodeFamiliesArr; -%newobject ParaMEDMEM::MEDFileUMesh::getMeshAtRank; -%newobject ParaMEDMEM::MEDFileUMesh::getRank0Mesh; -%newobject ParaMEDMEM::MEDFileUMesh::getRankM1Mesh; -%newobject ParaMEDMEM::MEDFileUMesh::getRankM2Mesh; -%newobject ParaMEDMEM::MEDFileUMesh::getRankM3Mesh; +%newobject ParaMEDMEM::MEDFileUMesh::getMeshAtLevel; +%newobject ParaMEDMEM::MEDFileUMesh::getLevel0Mesh; +%newobject ParaMEDMEM::MEDFileUMesh::getLevelM1Mesh; +%newobject ParaMEDMEM::MEDFileUMesh::getLevelM2Mesh; +%newobject ParaMEDMEM::MEDFileUMesh::getLevelM3Mesh; class MEDLoader { @@ -158,19 +158,18 @@ public: } static void WriteUMeshesPartition(const char *fileName, const char *meshName, PyObject *li, bool writeFromScratch) throw(INTERP_KERNEL::Exception) { - std::vector v=convertFieldDoubleVecFromPy(li); + std::vector v=convertUMeshVecFromPy(li); MEDLoader::WriteUMeshesPartition(fileName,meshName,v,writeFromScratch); } static void WriteUMeshesPartitionDep(const char *fileName, const char *meshName, PyObject *li, bool writeFromScratch) throw(INTERP_KERNEL::Exception) { - std::vector v=convertFieldDoubleVecFromPy(li); + std::vector v=convertUMeshVecFromPy(li); MEDLoader::WriteUMeshesPartitionDep(fileName,meshName,v,writeFromScratch); } static void WriteUMeshes(const char *fileName, PyObject *li, bool writeFromScratch) throw(INTERP_KERNEL::Exception) { - std::vector v=convertFieldDoubleVecFromPy(li); - std::vector v2(v.begin(),v.end()); - MEDLoader::WriteUMeshes(fileName,v2,writeFromScratch); + std::vector v=convertUMeshVecFromPy(li); + MEDLoader::WriteUMeshes(fileName,v,writeFromScratch); } static PyObject *GetTypesOfField(const char *fileName, const char *fieldName, const char *meshName) throw(INTERP_KERNEL::Exception) { @@ -210,3 +209,12 @@ public: %include "MEDFileMesh.hxx" +%extend ParaMEDMEM::MEDFileUMesh +{ + void setGroupsAtLevel(int meshDimRelToMaxExt, PyObject *li, bool renum=true) throw(INTERP_KERNEL::Exception) + { + std::vector grps; + convertPyObjToVecDataArrayIntCst(li,grps); + self->setGroupsAtLevel(meshDimRelToMaxExt,grps,renum); + } +} diff --git a/src/MEDLoader/Swig/MEDLoaderTest3.py b/src/MEDLoader/Swig/MEDLoaderTest3.py index 57fa71803..252a668ad 100644 --- a/src/MEDLoader/Swig/MEDLoaderTest3.py +++ b/src/MEDLoader/Swig/MEDLoaderTest3.py @@ -29,10 +29,10 @@ class MEDLoaderTest(unittest.TestCase): mname="ExampleOfMultiDimW" medmesh=MEDFileUMesh.New(fileName,mname) self.assertEqual((0,-1),medmesh.getNonEmptyLevels()) - m1_0=medmesh.getRank0Mesh() + m1_0=medmesh.getLevel0Mesh() m1_1=MEDLoader.ReadUMeshFromFile(fileName,mname,0) self.assertTrue(m1_0.isEqual(m1_1,1e-12)); - m2_0=medmesh.getRankM1Mesh() + m2_0=medmesh.getLevelM1Mesh() m2_1=MEDLoader.ReadUMeshFromFile(fileName,mname,-1) self.assertTrue(m2_0.isEqual(m2_1,1e-12)); pass @@ -42,7 +42,7 @@ class MEDLoaderTest(unittest.TestCase): outFileName="MEDFileMesh1.med" medmesh=MEDFileUMesh.New(fileName,mname) self.assertEqual((0,),medmesh.getNonEmptyLevels()) - m1_0=medmesh.getRank0Mesh() + m1_0=medmesh.getLevel0Mesh() m1_1=MEDLoader.ReadUMeshFromFile(fileName,mname,0) self.assertTrue(m1_0.isEqual(m1_1,1e-12)); g1_0=medmesh.getGroup(0,"mesh2") @@ -70,21 +70,179 @@ class MEDLoaderTest(unittest.TestCase): self.assertEqual([19,2,3,4,5,14,15,16],medmesh.getGroupsArr(0,["mesh2","mesh4","mesh3"]).getValues()); famn=medmesh.getFamilyNameGivenId(0) self.assertEqual(range(60),medmesh.getNodeFamilyArr(famn).getValues()); - # without renum + #without renum self.assertEqual([2,3,5,14,16],medmesh.getGroupArr(0,"mesh2",False).getValues()); self.assertEqual([2,3,16],medmesh.getFamilyArr(0,"Family_2",False).getValues()); self.assertEqual([2,3,5,14,16],medmesh.getFamiliesArr(0,["Family_4","Family_2"],False).getValues()); self.assertEqual([0,2,3,4,5,14,15,16],medmesh.getGroupsArr(0,["mesh2","mesh3","mesh4"],False).getValues()); self.assertEqual(range(60),medmesh.getNodeFamilyArr(famn,False).getValues()); pass + + # this tests emulates MEDMEM ( Except that it works ! ) The permutation are NOT taken into account def testMEDMesh3(self): + outFileName="MEDFileMesh3.med" + c=DataArrayDouble.New() + coords=[-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 ]; + targetConn=[0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4] + c.setValues(coords,9,2) + m=MEDCouplingUMesh.New(); + m.setMeshDimension(2); + m.allocateCells(5); + m.insertNextCell(NORM_TRI3,3,targetConn[4:7]) + m.insertNextCell(NORM_TRI3,3,targetConn[7:10]) + m.insertNextCell(NORM_QUAD4,4,targetConn[0:4]) + m.insertNextCell(NORM_QUAD4,4,targetConn[10:14]) + m.insertNextCell(NORM_QUAD4,4,targetConn[14:18]) + m.finishInsertingCells(); + m.setCoords(c) + m.checkCoherency() + m1=MEDCouplingUMesh.New(); + m1.setMeshDimension(1); + m1.allocateCells(3); + m1.insertNextCell(NORM_SEG2,2,[1,4]) + m1.insertNextCell(NORM_SEG2,2,[3,6]) + m1.insertNextCell(NORM_SEG3,3,[2,8,5]) + m1.finishInsertingCells(); + m1.setCoords(c) + m1.checkCoherency() + m2=MEDCouplingUMesh.New(); + m2.setMeshDimension(0); + m2.allocateCells(4); + m2.insertNextCell(NORM_POINT1,1,[1]) + m2.insertNextCell(NORM_POINT1,1,[3]) + m2.insertNextCell(NORM_POINT1,1,[2]) + m2.insertNextCell(NORM_POINT1,1,[6]) + m2.finishInsertingCells(); + m2.setCoords(c) + m2.checkCoherency() + # mm=MEDFileUMesh.New() mm.setName("MyFirstMEDCouplingMEDmesh") mm.setDescription("IHopeToConvinceLastMEDMEMUsers") + mm.setCoords(c) + mm.setMeshAtLevelOld(-1,m1); + mm.setMeshAtLevelOld(0,m); + mm.setMeshAtLevelOld(-2,m2); + # playing with groups + g1_2=DataArrayInt.New() + g1_2.setValues([1,3],2,1) + g1_2.setName("G1") + g2_2=DataArrayInt.New() + g2_2.setValues([1,2,3],3,1) + g2_2.setName("G2") + mm.setGroupsAtLevel(0,[g1_2,g2_2],False) + g1_1=DataArrayInt.New() + g1_1.setValues([0,1,2],3,1) + g1_1.setName("G1") + g2_1=DataArrayInt.New() + g2_1.setValues([0,2],2,1) + g2_1.setName("G2") + mm.setGroupsAtLevel(-1,[g1_1,g2_1],False) + g1_N=DataArrayInt.New() + g1_N.setValues(range(8),8,1) + g1_N.setName("G1") + g2_N=DataArrayInt.New() + g2_N.setValues(range(9),9,1) + g2_N.setName("G2") + mm.setGroupsAtLevel(1,[g1_N,g2_N],False) + # check content of mm + t=mm.getGroupArr(0,"G1",False) + self.assertTrue(g1_2.isEqual(t)); + t=mm.getGroupArr(0,"G2",False) + self.assertTrue(g2_2.isEqual(t)); + t=mm.getGroupArr(-1,"G1",False) + self.assertTrue(g1_1.isEqual(t)); + t=mm.getGroupArr(-1,"G2",False) + self.assertTrue(g2_1.isEqual(t)); + t=mm.getGroupArr(1,"G1",False) + self.assertTrue(g1_N.isEqual(t)); + t=mm.getGroupArr(1,"G2",False) + self.assertTrue(g2_N.isEqual(t)); + # + mm.write(outFileName,2); + pass + + # this test is the testMEDMesh3 except that permutation is dealed here + def testMEDMesh4(self): + outFileName="MEDFileMesh4.med" c=DataArrayDouble.New() coords=[-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 ]; + targetConn=[0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4] c.setValues(coords,9,2) + m=MEDCouplingUMesh.New(); + m.setMeshDimension(2); + m.allocateCells(5); + m.insertNextCell(NORM_QUAD4,4,targetConn[0:4]) + m.insertNextCell(NORM_TRI3,3,targetConn[4:7]) + m.insertNextCell(NORM_TRI3,3,targetConn[7:10]) + m.insertNextCell(NORM_QUAD4,4,targetConn[10:14]) + m.insertNextCell(NORM_QUAD4,4,targetConn[14:18]) + m.finishInsertingCells(); + m.setCoords(c) + m.checkCoherency() + m1=MEDCouplingUMesh.New(); + m1.setMeshDimension(1); + m1.allocateCells(3); + m1.insertNextCell(NORM_SEG2,2,[1,4]) + m1.insertNextCell(NORM_SEG3,3,[2,8,5]) + m1.insertNextCell(NORM_SEG2,2,[3,6]) + m1.finishInsertingCells(); + m1.setCoords(c) + m1.checkCoherency() + m2=MEDCouplingUMesh.New(); + m2.setMeshDimension(0); + m2.allocateCells(4); + m2.insertNextCell(NORM_POINT1,1,[1]) + m2.insertNextCell(NORM_POINT1,1,[3]) + m2.insertNextCell(NORM_POINT1,1,[2]) + m2.insertNextCell(NORM_POINT1,1,[6]) + m2.finishInsertingCells(); + m2.setCoords(c) + m2.checkCoherency() + # + mm=MEDFileUMesh.New() + mm.setName("My2ndMEDCouplingMEDmesh") + mm.setDescription("ThisIsImpossibleToDoWithMEDMEM") mm.setCoords(c) + renumNode=DataArrayInt.New() + renumNode.setValues([10,11,12,13,14,15,16,17,18],9,1) + mm.setRenumArr(1,renumNode) + mm.setMeshAtLevel(-1,m1); + mm.setMeshAtLevel(0,m); + mm.setMeshAtLevel(-2,m2); + # playing with groups + g1_2=DataArrayInt.New() + g1_2.setValues([2,3],2,1) + g1_2.setName("G1") + g2_2=DataArrayInt.New() + g2_2.setValues([2,0,3],3,1) + g2_2.setName("G2") + mm.setGroupsAtLevel(0,[g1_2,g2_2],True) + g1_1=DataArrayInt.New() + g1_1.setValues([0,2,1],3,1) + g1_1.setName("G1") + g2_1=DataArrayInt.New() + g2_1.setValues([0,2],2,1) + g2_1.setName("G2") + mm.setGroupsAtLevel(-1,[g1_1,g2_1],True) + g1_N=DataArrayInt.New() + g1_N.setValues([10,11,12,13,14,15,16,17],8,1) + g1_N.setName("G1") + g2_N=DataArrayInt.New() + g2_N.setValues([10,11,12,13,14,15,16,17,18],9,1) + g2_N.setName("G2") + mm.setGroupsAtLevel(1,[g1_N,g2_N],True) + # check content of mm + t=mm.getGroupArr(0,"G1",True) + self.assertTrue(g1_2.isEqual(t)); + t=mm.getGroupArr(0,"G2",True) + self.assertTrue(g2_2.isEqual(t)); + t=mm.getGroupArr(-1,"G1",True) + self.assertTrue(g1_1.isEqual(t)); + t=mm.getGroupArr(-1,"G2",True) + self.assertTrue(g2_1.isEqual(t)); + # + mm.write(outFileName,2); pass pass diff --git a/src/MEDPartitioner/MEDPARTITIONER_MESHCollection.cxx b/src/MEDPartitioner/MEDPARTITIONER_MESHCollection.cxx index 133e5f951..3865e8e7c 100644 --- a/src/MEDPartitioner/MEDPARTITIONER_MESHCollection.cxx +++ b/src/MEDPartitioner/MEDPARTITIONER_MESHCollection.cxx @@ -23,7 +23,7 @@ #include "MEDCouplingNormalizedUnstructuredMesh.hxx" #include "MEDCouplingMemArray.hxx" #include "PointLocator3DIntersectorP0P0.hxx" - +#include "BBTree.txx" #include "MEDPARTITIONER_utils.hxx" #include "MEDPARTITIONER_Graph.hxx" @@ -39,7 +39,7 @@ #include "MEDPARTITIONER_MESHCollectionDriver.hxx" #include "MEDPARTITIONER_MESHCollectionMedXMLDriver.hxx" #include "MEDPARTITIONER_MESHCollectionMedAsciiDriver.hxx" - +#include "MEDPARTITIONER_JointFinder.hxx" #include "MEDPARTITIONER_UserGraph.hxx" @@ -52,6 +52,7 @@ #include #include +#include #ifndef WNT # include @@ -111,66 +112,10 @@ MESHCollection::MESHCollection(MESHCollection& initialCollection, Topology* topo _create_empty_groups(create_empty_groups) { - _mesh.resize(_topology->nbDomain()); - - //splitting the initial domains into smaller bits - - std::vector > splitMeshes; - splitMeshes.resize(topology->nbDomain()); - for (int inew=0; inewnbDomain();inew++) - splitMeshes[inew].resize(initialCollection.getTopology()->nbDomain()); - std::vector > > new2oldIds(initialCollection.getTopology()->nbDomain()); + castCellMeshes(initialCollection, new2oldIds); - for (int iold=0; ioldnbDomain();iold++) - { - if (!isParallelMode() || initialCollection._domain_selector->isMyDomain(iold)) - { - int size=(initialCollection._mesh)[iold]->getNumberOfCells(); - std::vector globalids(size); - initialCollection.getTopology()->getCellList(iold, &globalids[0]); - std::vector ilocalnew(size); - std::vector ipnew(size); - topology->convertGlobalCellList(&globalids[0],size,&ilocalnew[0],&ipnew[0]); - new2oldIds[iold].resize(topology->nbDomain()); - for (int i=0; inbDomain();inew++) - { - // cout <<"inew:"<buildPartOfMySelf(&new2oldIds[iold][inew][0],&new2oldIds[iold][inew][0]+new2oldIds[iold][inew].size(),true); - // cout <<"small domain "<getNumberOfCells()<<" cells"<nbDomain() ;inew++) - { - std::vector meshes(splitMeshes[inew].size()); - for (int i=0; i< splitMeshes[inew].size();i++) - meshes[i]=splitMeshes[inew][i]; - - if (!isParallelMode()||_domain_selector->isMyDomain(inew)) - { - _mesh[inew]=ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes(meshes); - _mesh[inew]->zipCoords(); - - // std::cout<<"creating mesh #"<getNumberOfCells()<<" cells and "<<_mesh[inew]->getNumberOfNodes()<<" nodes"<decrRef(); - } //casting cell families on new meshes _cellFamilyIds.resize(topology->nbDomain()); castIntField(initialCollection.getMesh(), this->getMesh(),initialCollection.getCellFamilyIds(),_cellFamilyIds, new2oldIds); @@ -182,22 +127,30 @@ MESHCollection::MESHCollection(MESHCollection& initialCollection, Topology* topo // treating faces ///////////////// - std::multimap, pair > nodeMapping; + NodeMapping nodeMapping; createNodeMapping(initialCollection, nodeMapping); std::vector > > new2oldFaceIds; castMeshes(initialCollection.getFaceMesh(), this->getFaceMesh(),initialCollection, nodeMapping, new2oldFaceIds); + + //////////////////// + //treating families + //////////////////// + _faceFamilyIds.resize(topology->nbDomain()); //allocating family ids arrays for (int inew=0; inewnbDomain();inew++) { + if(isParallelMode() && !_domain_selector->isMyDomain(inew)) continue; _cellFamilyIds[inew]=ParaMEDMEM::DataArrayInt::New(); int* ptrCellIds=new int[_mesh[inew]->getNumberOfCells()]; + for (int i=0; i< _mesh[inew]->getNumberOfCells();i++) ptrCellIds[i]=0; _cellFamilyIds[inew]->useArray(ptrCellIds,true, ParaMEDMEM::CPP_DEALLOC,_mesh[inew]->getNumberOfCells(),1); _faceFamilyIds[inew]=ParaMEDMEM::DataArrayInt::New(); int* ptrFaceIds=new int[_faceMesh[inew]->getNumberOfCells()]; + for (int i=0; i< _faceMesh[inew]->getNumberOfCells();i++) ptrFaceIds[i]=0; _faceFamilyIds[inew]->useArray(ptrFaceIds,true, ParaMEDMEM::CPP_DEALLOC,_faceMesh[inew]->getNumberOfCells(),1); } @@ -212,36 +165,159 @@ MESHCollection::MESHCollection(MESHCollection& initialCollection, Topology* topo } +/*! +Creates the meshes using the topology underlying he mesh collection and the mesh data coming from the ancient collection +\param initialCollection collection from which the data is extracted to create the new meshes + */ + +void MESHCollection::castCellMeshes(MESHCollection& initialCollection, std::vector > >& new2oldIds) +{ + if (_topology==0) throw INTERP_KERNEL::Exception("Topology has not been defined on call to castCellMeshes"); + _mesh.resize(_topology->nbDomain()); + + //splitting the initial domains into smaller bits + + std::vector > splitMeshes; + splitMeshes.resize(_topology->nbDomain()); + for (int inew=0; inew<_topology->nbDomain();inew++) + { + splitMeshes[inew].resize(initialCollection.getTopology()->nbDomain()); + std::fill(&(splitMeshes[inew][0]),&(splitMeshes[inew][0])+splitMeshes[inew].size(),(ParaMEDMEM::MEDCouplingUMesh*)0); + } + + for (int iold=0; ioldnbDomain();iold++) + { + if (!isParallelMode() || initialCollection._domain_selector->isMyDomain(iold)) + { + int size=(initialCollection._mesh)[iold]->getNumberOfCells(); + std::vector globalids(size); + initialCollection.getTopology()->getCellList(iold, &globalids[0]); + std::vector ilocalnew(size); + std::vector ipnew(size); + _topology->convertGlobalCellList(&globalids[0],size,&ilocalnew[0],&ipnew[0]); + new2oldIds[iold].resize(_topology->nbDomain()); + for (int i=0; inbDomain();inew++) + { + splitMeshes[inew][iold]=(ParaMEDMEM::MEDCouplingUMesh*)(initialCollection.getMesh())[iold]->buildPartOfMySelf(&new2oldIds[iold][inew][0],&new2oldIds[iold][inew][0]+new2oldIds[iold][inew].size(),true); + } + } + } + + if (isParallelMode()) + { + for (int iold=0; ioldnbDomain();iold++) + for(int inew=0;inew<_topology->nbDomain();inew++) + { + if (initialCollection._domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)) continue; + if(initialCollection._domain_selector->isMyDomain(iold)) + { + _domain_selector->sendMesh(*(splitMeshes[inew][iold]),_domain_selector->getProcessorID(inew)); + std::cout<<"send iold"<getNumberOfCells()<isMyDomain(inew)) + { + _domain_selector->recvMesh(splitMeshes[inew][iold],_domain_selector->getProcessorID(iold)); + std::cout<<"recv iold"<nbDomain() ;inew++) + { + std::vector meshes; + + for (int i=0; i< splitMeshes[inew].size();i++) + if (splitMeshes[inew][i]!=0) meshes.push_back(splitMeshes[inew][i]); + std::cout<< "nb of meshes"<isMyDomain(inew)) + { + _mesh[inew]=ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes(meshes); + _mesh[inew]->zipCoords(); + } + for (int i=0; i< splitMeshes[inew].size();i++) + if (splitMeshes[inew][i]!=0) splitMeshes[inew][i]->decrRef(); + } + +} + /*! \param initialCollection source mesh collection \param nodeMapping structure containing the correspondency between nodes in the initial collection and the node(s) in the new collection */ -void MESHCollection::createNodeMapping( MESHCollection& initialCollection, std::multimap,pair >& nodeMapping) +void MESHCollection::createNodeMapping( MESHCollection& initialCollection, NodeMapping& nodeMapping) { + + // NodeMapping reverseNodeMapping; for (int iold=0; ioldnbDomain();iold++) { - std::map >, int > nodeClassifier; - for (int inode=0; inodegetNumberOfNodes(); inode++) + + double* bbox; + BBTree<3>* tree; + if (!isParallelMode() || (_domain_selector->isMyDomain(iold))) { + // std::map >, int > nodeClassifier; + int nvertices=getMesh(iold)->getNumberOfNodes(); + bbox=new double[nvertices*6]; ParaMEDMEM::DataArrayDouble* coords = initialCollection.getMesh(iold)->getCoords(); - double* coordsPtr=coords->getPointer()+initialCollection.getMesh(iold)->getMeshDimension()*inode; - pair > tripleDouble =std::make_pair(coordsPtr[0],std::make_pair(coordsPtr[1],coordsPtr[2])); - nodeClassifier.insert(make_pair(tripleDouble,inode)); - + double* coordsPtr=coords->getPointer(); + for (int i=0; i(bbox,0,0,nvertices,1e-12); } + for (int inew=0; inew<_topology->nbDomain(); inew++) { - for (int inode=0; inode<_mesh[inew]->getNumberOfNodes();inode++) + //sending meshes for parallel computation + if (isParallelMode() && _domain_selector->isMyDomain(inew) && !_domain_selector->isMyDomain(iold)) { - ParaMEDMEM::DataArrayDouble* coords = getMesh(inew)->getCoords(); - double* coordsPtr=coords->getPointer()+inode*getMesh(inew)->getMeshDimension(); - std::map >, int >::iterator iter= - nodeClassifier.find(make_pair(coordsPtr[0],make_pair(coordsPtr[1],coordsPtr[2]))); - nodeMapping.insert(make_pair(make_pair(iold,iter->second),make_pair(inew,inode))); - // cout <<"("<second<<")-->("<getProcessorID(iold)<sendMesh(*(getMesh(inew)), _domain_selector->getProcessorID(iold)); + + } + else if (isParallelMode() && !_domain_selector->isMyDomain(inew)&& _domain_selector->isMyDomain(iold)) + { + ParaMEDMEM::MEDCouplingUMesh* mesh; + std::cout<<"recvFrom" << _domain_selector->getProcessorID(inew) <recvMesh(mesh, _domain_selector->getProcessorID(inew)); + ParaMEDMEM::DataArrayDouble* coords = mesh->getCoords(); + for (int inode=0; inodegetNumberOfNodes();inode++) + { + double* coordsPtr=coords->getPointer()+inode*3; + std::vector elems; + tree->getElementsAroundPoint(coordsPtr,elems); + if (elems.size()==0) continue; + nodeMapping.insert(make_pair(make_pair(iold,elems[0]),make_pair(inew,inode))); + } + } + else if (!isParallelMode() || (_domain_selector->isMyDomain(inew) && _domain_selector->isMyDomain(iold))) + { + ParaMEDMEM::DataArrayDouble* coords = getMesh(inew)->getCoords(); + + for (int inode=0; inode<_mesh[inew]->getNumberOfNodes();inode++) + { + + double* coordsPtr=coords->getPointer()+inode*3; + std::vector elems; + tree->getElementsAroundPoint(coordsPtr,elems); + if (elems.size()==0) {continue;} + nodeMapping.insert(make_pair(make_pair(iold,elems[0]),make_pair(inew,inode))); + cout << "inode :" <("<& meshesCastFrom,std::vector& meshesCastTo, MESHCollection& initialCollection,const multimap,pair >& nodeMapping, std::vector > >& new2oldIds) +void MESHCollection::castMeshes(std::vector& meshesCastFrom,std::vector& meshesCastTo, MESHCollection& initialCollection,const NodeMapping& nodeMapping, std::vector > >& new2oldIds) { //splitMeshes structure will contain the partition of @@ -267,11 +343,13 @@ void MESHCollection::castMeshes(std::vector& mesh new2oldIds.resize(meshesCastFrom.size()); + //loop over the old domains to analyse the faces and decide //on which new domain they belong for (int iold=0; ioldisMyDomain(iold)) continue; new2oldIds[iold].resize(newSize); for (int ielem=0;ielemgetNumberOfCells();ielem++) { @@ -284,7 +362,9 @@ void MESHCollection::castMeshes(std::vector& mesh for (int inode=0;inode,pair >::const_iterator MI; - pair myRange = nodeMapping.equal_range(make_pair(iold,nodes[inode])); + int mynode=nodes[inode]; + if (mynode <0 || mynode > 1000000000) exit(1); + pair myRange = nodeMapping.equal_range(make_pair(iold,mynode)); // cout << iold <<" " <& mesh //creating the splitMeshes from the face ids for (int inew=0;inew<_topology->nbDomain();inew++) { - cout<<"nb faces "<buildPartOfMySelf(&new2oldIds[iold][inew][0],&new2oldIds[iold][inew][0]+new2oldIds[iold][inew].size(),true)); } } // send/receive stuff + if (isParallelMode()) + for (int iold=0; ioldisMyDomain(iold) && !_domain_selector->isMyDomain(inew)) + _domain_selector->sendMesh(*(splitMeshes[inew][iold]), _domain_selector->getProcessorID(inew)); + if (!_domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)) + _domain_selector->recvMesh(splitMeshes[inew][iold], _domain_selector->getProcessorID(iold)); + } //recollecting the bits of splitMeshes to fuse them into one meshesCastTo.resize(newSize); for (int inew=0; inew < newSize;inew++) { - vector myMeshes(meshesCastFrom.size()); + vector myMeshes; for (int iold=0; iold < meshesCastFrom.size();iold++) { - myMeshes[iold]=splitMeshes[inew][iold]; - cout<<"number of nodes"<getNumberOfNodes()<getNumberOfCells()<getMeshDimension()<checkCoherency(); meshesCastTo[inew]->zipCoords(); - meshesCastTo[inew]->checkCoherency(); for (int iold=0; iold < meshesCastFrom.size();iold++) - splitMeshes[inew][iold]->decrRef(); - cout<<"creating face mesh #"<getNumberOfCells()<<" cells and "<getNumberOfNodes()<<" nodes"<decrRef(); } } @@ -346,6 +429,7 @@ void MESHCollection::castIntField(std::vector& me vector splitIds; for (int iold=0; iold < meshesCastFrom.size();iold++) { + if (isParallelMode() && !_domain_selector->isMyDomain(iold)) continue; int* ptr=&(new2oldMapping[iold][inew][0]); splitIds.push_back(arrayFrom[iold]->selectByTupleId(ptr,ptr+new2oldMapping[iold][inew].size())); } @@ -561,7 +645,7 @@ vector& MESHCollection::getFaceMesh() { return _faceMesh; } -ParaMEDMEM::MEDCouplingUMesh* MESHCollection::getMesh(int idomain) +ParaMEDMEM::MEDCouplingUMesh* MESHCollection::getMesh(int idomain) const { return _mesh[idomain]; } @@ -604,9 +688,25 @@ void MESHCollection::buildCellGraph(MEDPARTITIONER::MEDSKYLINEARRAY* & array,int multimap< int, int > node2cell; multimap< int, int > cell2cell; + vector > > commonDistantNodes; + int nbdomain=_topology->nbDomain(); + if (isParallelMode()) + { + _joint_finder=new JointFinder(*this); + _joint_finder->findCommonDistantNodes(); + commonDistantNodes=_joint_finder->getDistantNodeCell(); + } + //looking for reverse nodal connectivity i global numbering - for (int idomain=0; idomain<_topology->nbDomain();idomain++) + for (int idomain=0; idomainisMyDomain(idomain)) continue; + int offset=0, procOffset=0; + if (isParallelMode()) + { + offset=_domain_selector->getDomainShift(idomain); + procOffset=_domain_selector->getProcShift()+1; + } ParaMEDMEM::DataArrayInt* index=ParaMEDMEM::DataArrayInt::New(); ParaMEDMEM::DataArrayInt* revConn=ParaMEDMEM::DataArrayInt::New(); _mesh[idomain]->getReverseNodalConnectivity(revConn,index); @@ -617,7 +717,7 @@ void MESHCollection::buildCellGraph(MEDPARTITIONER::MEDSKYLINEARRAY* & array,int index_ptr[_mesh[idomain]->getNumberOfNodes()], globalRevConn->getPointer()); - int min=10000000,max=-1; + int min=std::numeric_limits::max(),max=-1; for (int i=0; i< index_ptr[_mesh[idomain]->getNumberOfNodes()];i++) { if (revConn->getPointer()[i]getPointer()[i]; @@ -626,17 +726,38 @@ void MESHCollection::buildCellGraph(MEDPARTITIONER::MEDSKYLINEARRAY* & array,int cout <<"min:"<getNumberOfNodes()]; _topology->getNodeList(idomain,globalNodeIds); - + int* globalRevConnPtr=globalRevConn->getPointer(); for (int i=0; i<_mesh[idomain]->getNumberOfNodes();i++) { for (int icell=index_ptr[i]; icellnbDomain();mydomain++) + { + if (_domain_selector->isMyDomain(mydomain)) continue; + // for (int idomain=0; idomain<_topology->nbDomain();idomain++) + // { + // if (_domain_selector->isMyDomain(idomain)) continue; + multimap::iterator iter; + for (iter=commonDistantNodes[idomain][mydomain].begin();iter!=commonDistantNodes[idomain][mydomain].end();iter++) + { + int ilocnode=iter->first; + int icell=iter->second; + + node2cell.insert(make_pair(globalNodeIds[ilocnode],icell+offset)); + // cout<<"pair "<decrRef(); revConn->decrRef(); index->decrRef(); - delete[]globalNodeIds; + delete[] globalNodeIds; } //creating graph arcs (cell to cell relations) @@ -646,6 +767,21 @@ void MESHCollection::buildCellGraph(MEDPARTITIONER::MEDSKYLINEARRAY* & array,int // means 6 arcs (0,1), (0,2), (0,3), (1,2), (1,3), (2,3) // in present version arcs are not doubled but reflexive (1,1) arcs are present for each cell + int mincell,maxcell; + if (isParallelMode()) + { + mincell=_domain_selector->getProcShift()+1; + maxcell=mincell; + for (int i=0; iisMyDomain(i)) maxcell+=_mesh[i]->getNumberOfCells(); + + } + else + { + mincell=0; + maxcell=_topology->nbCells(); + } + cout<<"mincell"<nbNodes();inode++) { typedef multimap::const_iterator MI; @@ -654,7 +790,7 @@ void MESHCollection::buildCellGraph(MEDPARTITIONER::MEDSKYLINEARRAY* & array,int { for (MI cell2 = myRange.first; cell2!=myRange.second; cell2++) { - if (cell1->second!=cell2->second) cell2cell.insert(make_pair(cell1->second,cell2->second)); + if (cell1->second!=cell2->second&&cell1->second>=mincell&&cell1->secondsecond,cell2->second)); } } } @@ -682,204 +818,102 @@ void MESHCollection::buildCellGraph(MEDPARTITIONER::MEDSKYLINEARRAY* & array,int } array=new MEDPARTITIONER::MEDSKYLINEARRAY(index,value); - cout<<"taille du graphe créé "<getNumberOf()<nbDomain(); i++) - // { - // cell_number+=_topology->getCellNumber(i); - // node_number+=_topology->getNodeNumber(i); - // } - // //list of cells for a given node - // //vector< vector > node2cell(node_number); - // map< int, vector > node2cell; - - // //list of nodes for a given cell - // //vector< vector > cell2node(cell_number); - // map< int, vector > cell2node; - - // // map type_cell_list; - - // //tagging for the indivisible regions - // int* indivisible_tag=0; - // bool has_indivisible_regions=false; - // // if (!_indivisible_regions.empty()) - // // { - // // has_indivisible_regions=true; - // // indivisible_tag=new int[_topology->nbCells()]; - // // treatIndivisibleRegions(indivisible_tag); - // // } - - // fillGlobalConnectivity(node2cell, cell2node ); - - // cout << "beginning of skyline creation"<nbCells()]; - // int** temp=new int*[_topology->nbCells()]; - // int** temp_edgeweight=0; - // // if (has_indivisible_regions) - // // temp_edgeweight=new int*[_topology->nbCells()]; - - // int cell_glob_shift = 0; - - // // Get connection to cells on other procs - // multimap< int, int > loc2dist; // global cell ids on this proc -> other proc cells - // if ( isParallelMode() ) - // { - // cell_glob_shift = _domain_selector->getProcShift(); - - // set loc_domains; // domains on this proc - // for ( int idom = 0; idom < _mesh.size(); ++idom ) - // if ( _mesh[ idom ] ) - // loc_domains.insert( idom ); - - // for ( int idom = 0; idom < _mesh.size(); ++idom ) - // { - // if ( !_mesh[idom] ) continue; - // vector loc2glob_corr; // pairs of corresponding cells (loc_loc & glob_dist) - // retrieveDriver()->readLoc2GlobCellConnect(idom, loc_domains, _domain_selector, loc2glob_corr); - // //MEDMEM::STRING out; - // for ( int i = 0; i < loc2glob_corr.size(); i += 2 ) - // { - // int glob_here = _topology->convertCellToGlobal(idom,loc2glob_corr[i]); - // int glob_there = loc2glob_corr[i+1]; - // loc2dist.insert ( make_pair( glob_here, glob_there)); - // //out << glob_here << "-" << glob_there << " "; - // } - // //cout << "\nRank " << _domain_selector->rank() << ": BndCZ: " << out << endl; - // } - // } - - // //going across all cells - - // map cells_neighbours; - // for (int i=0; i< _topology->nbCells(); i++) - // { - - - // vector cells(50); - - // // /*// cout << "size cell2node "<::const_iterator iternode=cell2node[i+1].begin(); - // // iternode!=cell2node[i+1].end(); - // // iternode++) - // // { - // // int nodeid=*iternode; - // // // cout << "size node2cell "<::const_iterator iter=node2cell[nodeid].begin(); - // // iter != node2cell[nodeid].end(); - // // iter++) - // // cells_neighbours[*iter]++; - // // }*/ - - // for (int inode=0; inode< cell2node[i+1].size(); inode++) - // { - // int nodeid=cell2node[i+1][inode]; - // // cout << "size node2cell "<::const_iterator iter=cells_neighbours.begin(); iter != cells_neighbours.end(); iter++) - // { - // if (iter->second >= dimension && iter->first != i+1) - // { - // cells.push_back(iter->first + cell_glob_shift); - // // cells[isize++]=iter->first; - // } - // } - // // add neighbour cells from distant domains - // multimap< int, int >::iterator loc_dist = loc2dist.find( i+1 ); - // for (; loc_dist!=loc2dist.end() && loc_dist->first==( i+1 ); ++loc_dist ) - // cells.push_back( loc_dist->second ); - - // size[i]=cells.size(); - // // size[i]=isize; - - // //cout << cells.size()<::const_iterator iter=cells.begin(); iter!=cells.end();iter++) - // //for(int j=0; jfirst-1]; - // // int tag2 = indivisible_tag[*iter-1]; - // // if (tag1==tag2 && tag1!=0) - // // temp_edgeweight[i][itemp]=_topology->nbCells()*100000; - // // else - // // temp_edgeweight[i][itemp]=1; - // // } - // itemp++; - // } - // cells_neighbours.clear(); - // } - // cout <<"end of graph definition"<nbCells()+1]; - // index[0]=1; - // for (int i=0; i<_topology->nbCells(); i++) - // index[i+1]=index[i]+size[i]; - - // node2cell.clear(); - // cell2node.clear(); - // // if (indivisible_tag!=0) delete [] indivisible_tag; - - // //SKYLINEARRAY structure holding the cell graph - // array= new MEDMEM::MEDSKYLINEARRAY(_topology->nbCells(),index[_topology->nbCells()]-index[0]); - // array->setIndex(index); - - // for (int i=0; i<_topology->nbCells(); i++) - // { - // array->setI(i+1,temp[i]); - // delete[]temp[i]; - // } - - // {// DEBUG - // // MEDMEM::STRING out; - // // const int* index= array->getIndex(); - // // const int* value =array->getValue(); - // // out << "\nRank " << (_domain_selector?_domain_selector->rank():0) << ": Index: "; - // // for ( int i = 0; i <= array->getNumberOf(); ++i ) - // // out << index[i] << " "; - // // out << "\nRank " << (_domain_selector?_domain_selector->rank():0) << ": Value: "; - // // for ( int i = 0; i < array->getLength(); ++i ) - // // out << value[i] << " "; - // // cout << out << endl; - // } - // // if (has_indivisible_regions) - // // { - // // edgeweights=new int[array->getLength()]; - // // for (int i=0; i<_topology->nbCells(); i++) - // // { - // // for (int j=index[i]; j > >& commonDistantNodes) +// { +// int nbdomain=_topology->nbDomain(); +// commonDistantNodes.resize(nbdomain); +// for (int i=0; inbProcs(); +// vector* > bbtree(nbdomain); +// vector rev(nbdomain); +// vectorrevIndx(nbdomain); +// int meshDim; +// int spaceDim; + +// for (int mydomain=0;mydomainisMyDomain(mydomain)) continue; +// meshDim=_mesh[mydomain]->getMeshDimension(); +// spaceDim= _mesh[mydomain]->getSpaceDimension(); +// rev[mydomain] = ParaMEDMEM::DataArrayInt::New(); +// revIndx[mydomain] = ParaMEDMEM::DataArrayInt::New(); +// _mesh[mydomain]->getReverseNodalConnectivity(rev[mydomain],revIndx[mydomain]); +// double* bbx=new double[2*spaceDim*_mesh[mydomain]->getNumberOfNodes()]; +// for (int i=0; i<_mesh[mydomain]->getNumberOfNodes()*spaceDim;i++) +// { +// const double* coords=_mesh[mydomain]->getCoords()->getConstPointer(); +// bbx[2*i]=(coords[i])-1e-12; +// bbx[2*i+1]=bbx[2*i]+2e-12; +// } +// bbtree[mydomain]=new BBTree<3> (bbx,0,0,_mesh[mydomain]->getNumberOfNodes(),-1e-12); +// } +// for (int isource=0;isourceisMyDomain(isource)&&_domain_selector->isMyDomain(itarget)) continue; +// if (_domain_selector->isMyDomain(isource)) +// { +// //preparing data for treatment on target proc +// int targetProc = _domain_selector->getProccessorID(itarget); + +// std::vector vec(spaceDim*_mesh[isource]->getNumberOfNodes()); +// std::copy(_mesh[isource]->getCoords()->getConstPointer(),_mesh[isource]->getCoords()->getConstPointer()+_mesh[isource]->getNumberOfNodes()*spaceDim,&vec[0]); +// _domain_selector->sendDoubleVec (vec,targetProc); + +// //retrieving target data for storage in commonDistantNodes array +// vector localCorrespondency; +// _domain_selector->recvIntVec(localCorrespondency, targetProc); +// cout<<"size "<isMyDomain(itarget)) +// { +// //receiving data from source proc +// int sourceProc = isource%nbproc; +// std::vector recvVec; +// _domain_selector->recvDoubleVec(recvVec,sourceProc); +// std::map commonNodes; // (local nodes, distant nodes) list +// for (int inode=0; inode<(recvVec.size()/meshDim);inode++) +// { +// double* bbox=new double[2*spaceDim]; +// for (int i=0; i inodes; +// bbtree[itarget]->getIntersectingElems(bbox,inodes); +// delete[] bbox; + +// if (inodes.size()>0) commonNodes.insert(make_pair(inodes[0],inode)); +// } +// std::vector nodeCellCorrespondency; +// for (map::iterator iter=commonNodes.begin();iter!=commonNodes.end();iter++) +// { +// const int*revIndxPtr=revIndx[itarget]->getConstPointer(); +// const int*revPtr=rev[itarget]->getConstPointer(); +// for (int icell=revIndxPtr[iter->first];icellfirst+1];icell++) +// { +// nodeCellCorrespondency.push_back(iter->second); +// nodeCellCorrespondency.push_back(revPtr[icell]); +// } +// } +// _domain_selector->sendIntVec(nodeCellCorrespondency, sourceProc); +// } + +// } +// } + + /*! Creates the partition corresponding to the cell graph and the partition number * * \param nbdomain number of subdomains for the newly created graph @@ -1776,16 +1810,13 @@ Topology* MESHCollection::createPartition(const int* partition) // } // } -// void MESHCollection::castField(const MESHCollection& old_collection, const string& fieldname, int itnumber, int ordernumber) +// void MESHCollection::castFieldDouble(const MESHCollection& old_collection, const string& fieldname, int itnumber, int ordernumber) // { // int type=old_collection.getDriver()->getFieldType(fieldname); // char field_char[80]; // strcpy(field_char,fieldname.c_str()); - -// if (type ==0) -// castFields(old_collection, field_char, itnumber, ordernumber); -// else -// castFields(old_collection, field_char, itnumber, ordernumber); +// std::vector fieldTypes =MEDLoader::GetTypesOfField(fileName, fieldName,meshName); + // } // void MESHCollection::castAllFields(const MESHCollection& initial_collection) @@ -1794,18 +1825,18 @@ Topology* MESHCollection::createPartition(const int* partition) // vector iternumber; // vector ordernumber; // vector types; -// initial_collection.getDriver()->readFileStruct(field_names,iternumber,ordernumber,types); + +// string filename=initial_collection.getDriver()->getFilename(); +// field_names=MEDLoader::GetAllFieldNames(filename.c_str()); + +// readFileStruct(field_names,iternumber,ordernumber,types); // for (int i=0; i(initial_collection, field_char, iternumber[i], ordernumber[i]); -// else -// castFields(initial_collection, field_char, iternumber[i], ordernumber[i]); +// castFieldDouble(initial_collection, field_char, iternumber[i], ordernumber[i]); // } // } @@ -2211,6 +2242,7 @@ void MESHCollection::setDomainNames(const std::string& name) { ostringstream oss; oss<isMyDomain(i)) _mesh[i]->setName(oss.str().c_str()); } } diff --git a/src/MEDPartitioner/MEDPARTITIONER_MESHCollection.hxx b/src/MEDPartitioner/MEDPARTITIONER_MESHCollection.hxx index 766224c9b..72fb537f1 100644 --- a/src/MEDPartitioner/MEDPARTITIONER_MESHCollection.hxx +++ b/src/MEDPartitioner/MEDPARTITIONER_MESHCollection.hxx @@ -26,7 +26,7 @@ //#include "boost/shared_ptr.hpp" #include #include - +#include namespace ParaMEDMEM { class MEDCouplingUMesh; @@ -42,9 +42,11 @@ namespace MEDPARTITIONER class ParaDomainSelector; class MEDSKYLINEARRAY; class CONNECTZONE; - + class JointFinder; typedef enum{MedAscii, MedXML, Undefined} DriverType; + typedef std::multimap, std::pair > NodeMapping ; + typedef std::vector > NodeList; class MEDPARTITIONER_EXPORT MESHCollection { @@ -126,7 +128,7 @@ namespace MEDPARTITIONER std::vector& getFaceMesh() ; std::vector >&getGroupMeshes(); - ParaMEDMEM::MEDCouplingUMesh* getMesh(int idomain); + ParaMEDMEM::MEDCouplingUMesh* getMesh(int idomain) const; ParaMEDMEM::MEDCouplingUMesh* getFaceMesh(int idomain); std::vector& getGroupMeshes(int idomain); @@ -141,7 +143,7 @@ namespace MEDPARTITIONER //getting a pointer to topology Topology* getTopology() const ; - + ParaDomainSelector* getParaDomainSelector() const{return _domain_selector;} //settig a new topology void setTopology(Topology* topology); @@ -158,6 +160,7 @@ namespace MEDPARTITIONER //creates the node mapping between an old collection and the present one void createNodeMapping( MESHCollection& initialCollection, std::multimap,std::pair >& nodeMapping); + void castCellMeshes(MESHCollection& initialCollection, std::vector > >& new2oldIds); //creates faces on the new collection void castFaceMeshes( MESHCollection& initialCollection,const std::multimap,std::pair >& nodeMapping); @@ -212,6 +215,7 @@ namespace MEDPARTITIONER void castIntField(std::vector& meshesCastFrom,std::vector& meshesCastTo, std::vector& arrayFrom, std::vector& arrayTo, std::vector > > & old2newMapping); + void findCommonDistantNodes(std::vector > >& commonDistantNodes); // template // void castFields(const MESHCollection& old_collection, const string& fieldname, int itnumber, int ordernumber); @@ -278,6 +282,8 @@ namespace MEDPARTITIONER /*! flag specifying that groups must be created on all domains, even if they are empty*/ bool _create_empty_groups; + + JointFinder* _joint_finder; }; }//of namespace diff --git a/src/MEDPartitioner/MEDPARTITIONER_METISGraph.cxx b/src/MEDPartitioner/MEDPARTITIONER_METISGraph.cxx index 2c86aacb2..c7587119d 100644 --- a/src/MEDPartitioner/MEDPARTITIONER_METISGraph.cxx +++ b/src/MEDPartitioner/MEDPARTITIONER_METISGraph.cxx @@ -69,6 +69,7 @@ void METISGraph::partGraph(int ndomain, int edgecut; int* partition = new int[n]; + cout << "ParMETIS : n="<1) { if ( parallelizer ) @@ -77,7 +78,7 @@ void METISGraph::partGraph(int ndomain, // distribution of vertices of the graph among the processors int * vtxdist = parallelizer ? parallelizer->getNbVertOfProcs() : 0; MPI_Comm comm = MPI_COMM_WORLD; - + cout<<"vtxdist[1]"<<" "< ParaDomainSelector::gatherGraph(const Graph* graph) const // Make graph // ----------- - MEDMEM::MEDSKYLINEARRAY* array = - new MEDMEM::MEDSKYLINEARRAY( index_size-1, value_size, graph_index, graph_value, true ); +// MEDPARTITIONER::MEDSKYLINEARRAY* array = +// new MEDPARTITIONER::MEDSKYLINEARRAY( index_size-1, value_size, graph_index, graph_value, true ); - glob_graph = new UserGraph( array, partition, index_size-1 ); +// glob_graph = new UserGraph( array, partition, index_size-1 ); evaluateMemory(); @@ -535,7 +514,7 @@ void ParaDomainSelector::gatherNbCellPairs() // { // vector send_data, recv_data( joint->serialize( send_data )); -// int dest = getProccessorID( joint->distantDomain() ); +// int dest = getProcessorID( joint->distantDomain() ); // int tag = 1001 + jointId( joint->localDomain(), joint->distantDomain() ); // #ifdef HAVE_MPI2 @@ -583,7 +562,7 @@ int* ParaDomainSelector::exchangeSubentityIds( int loc_domain, int dist_domain, const vector& loc_ids_here ) const { int* loc_ids_dist = new int[ loc_ids_here.size()]; - int dest = getProccessorID( dist_domain ); + int dest = getProcessorID( dist_domain ); int tag = 2002 + jointId( loc_domain, dist_domain ); #ifdef HAVE_MPI2 MPI_Status status; @@ -621,7 +600,7 @@ int ParaDomainSelector::jointId( int local_domain, int distant_domain ) const // int ParaDomainSelector::getDomianOrder(int idomain, int nb_domains) const // { -// return nb_domains / nbProcs() * getProccessorID( idomain ) + idomain / nbProcs(); +// return nb_domains / nbProcs() * getProcessorID( idomain ) + idomain / nbProcs(); // } //================================================================================ @@ -665,3 +644,154 @@ int ParaDomainSelector::evaluateMemory() const } return _max_memory - _init_memory; } + + +void ParaDomainSelector::sendMesh(const ParaMEDMEM::MEDCouplingUMesh& mesh, int target) const +{ + + + // First stage : sending sizes + // ------------------------------ + vector tinyInfoLocal; + vector tinyInfoLocalS; + //Getting tiny info of local mesh to allow the distant proc to initialize and allocate + //the transmitted mesh. + mesh.getTinySerializationInformation(tinyInfoLocal,tinyInfoLocalS); + tinyInfoLocal.push_back(mesh.getNumberOfCells()); + #ifdef HAVE_MPI2 + int tinySize=tinyInfoLocal.size(); + MPI_Send(&tinySize, 1, MPI_INT, target, 1113, MPI_COMM_WORLD); + MPI_Send(&tinyInfoLocal[0], tinyInfoLocal.size(), MPI_INT, target, 1112, + MPI_COMM_WORLD); +#endif + ParaMEDMEM::DataArrayInt *v1Local=0; + ParaMEDMEM::DataArrayDouble *v2Local=0; + //serialization of local mesh to send data to distant proc. + mesh.serialize(v1Local,v2Local); + int nbLocalElems; + int* ptLocal=0; + if(v1Local) + { + nbLocalElems=v1Local->getNbOfElems(); + ptLocal=v1Local->getPointer(); + } + #ifdef HAVE_MPI2 + MPI_Send(ptLocal, nbLocalElems, MPI_INT, + target, 1111, MPI_COMM_WORLD); + #endif + double *ptLocal2=0; + if(v2Local) + { + nbLocalElems=v2Local->getNbOfElems(); + ptLocal2=v2Local->getPointer(); + } +#ifdef HAVE_MPI2 + MPI_Send(ptLocal2, nbLocalElems, MPI_DOUBLE, + target, 1110, MPI_COMM_WORLD); +#endif + if(v1Local) + v1Local->decrRef(); + if(v2Local) + v2Local->decrRef(); + +} +void ParaDomainSelector::recvMesh(ParaMEDMEM::MEDCouplingUMesh*& mesh, int source)const +{ + + // First stage : exchanging sizes + // ------------------------------ + vector tinyInfoDistant; + vector tinyInfoLocalS; + //Getting tiny info of local mesh to allow the distant proc to initialize and allocate + //the transmitted mesh. +#ifdef HAVE_MPI2 + MPI_Status status; + int tinyVecSize; + MPI_Recv(&tinyVecSize, 1, MPI_INT,source,1113,MPI_COMM_WORLD, &status); + tinyInfoDistant.resize(tinyVecSize); +#endif + std::fill(tinyInfoDistant.begin(),tinyInfoDistant.end(),0); + +#ifdef HAVE_MPI2 + MPI_Recv(&tinyInfoDistant[0], tinyVecSize, MPI_INT,source,1112,MPI_COMM_WORLD, &status); +#endif + ParaMEDMEM::DataArrayInt *v1Distant=ParaMEDMEM::DataArrayInt::New(); + ParaMEDMEM::DataArrayDouble *v2Distant=ParaMEDMEM::DataArrayDouble::New(); + //Building the right instance of copy of distant mesh. + ParaMEDMEM::MEDCouplingPointSet *distant_mesh_tmp=ParaMEDMEM::MEDCouplingPointSet::BuildInstanceFromMeshType((ParaMEDMEM::MEDCouplingMeshType)tinyInfoDistant[0]); + std::vector unusedTinyDistantSts; + mesh=dynamic_cast (distant_mesh_tmp); + + mesh->resizeForUnserialization(tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts); + int nbDistElem=0; + int *ptDist=0; + if(v1Distant) + { + nbDistElem=v1Distant->getNbOfElems(); + ptDist=v1Distant->getPointer(); + } +#ifdef HAVE_MPI2 + MPI_Recv(ptDist, nbDistElem, MPI_INT, + source,1111, + MPI_COMM_WORLD, &status); +#endif + int nbLocalElems=0; + double *ptDist2=0; + nbDistElem=0; + if(v2Distant) + { + nbDistElem=v2Distant->getNbOfElems(); + ptDist2=v2Distant->getPointer(); + } +#ifdef HAVE_MPI2 + MPI_Recv(ptDist2, nbDistElem, MPI_DOUBLE,source, 1110, MPI_COMM_WORLD, &status); +#endif + // + //mesh=dynamic_cast (distant_mesh_tmp); + //finish unserialization + mesh->unserialization(tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts); + std::cout<<"mesh size on recv"<getNumberOfCells()<decrRef(); + if(v2Distant) + v2Distant->decrRef(); + +} +void ParaDomainSelector::sendDoubleVec(const std::vector& vec, int target)const +{ + int size=vec.size(); +#ifdef HAVE_MPI2 + MPI_Send(&size,1,MPI_INT,target,1211, MPI_COMM_WORLD); + MPI_Send(const_cast(&vec[0]), size,MPI_DOUBLE, target, 1212, MPI_COMM_WORLD); +#endif +} +void ParaDomainSelector::recvDoubleVec(std::vector& vec, int source)const +{ + int size; +#ifdef HAVE_MPI2 + MPI_Status status; + MPI_Recv(&size,1,MPI_INT,source,1211, MPI_COMM_WORLD, &status); + vec.resize(size); + MPI_Recv(&vec[0],size,MPI_DOUBLE,source, 1212, MPI_COMM_WORLD,&status); +#endif +} +void ParaDomainSelector::sendIntVec(const std::vector& vec, int target)const +{ + int size=vec.size(); +#ifdef HAVE_MPI2 + MPI_Send(&size,1,MPI_INT,target,1211, MPI_COMM_WORLD); + MPI_Send(const_cast(&vec[0]), size,MPI_INT, target, 1212, MPI_COMM_WORLD); +#endif +} +void ParaDomainSelector::recvIntVec(std::vector& vec, int source)const +{ + int size; +#ifdef HAVE_MPI2 + MPI_Status status; + MPI_Recv(&size,1,MPI_INT,source,1211, MPI_COMM_WORLD, &status); + vec.resize(size); + MPI_Recv(&vec[0],size,MPI_INT,source, 1212, MPI_COMM_WORLD,&status); +#endif +} + diff --git a/src/MEDPartitioner/MEDPARTITIONER_ParaDomainSelector.hxx b/src/MEDPartitioner/MEDPARTITIONER_ParaDomainSelector.hxx index ee502b4b2..79c277502 100644 --- a/src/MEDPartitioner/MEDPARTITIONER_ParaDomainSelector.hxx +++ b/src/MEDPartitioner/MEDPARTITIONER_ParaDomainSelector.hxx @@ -65,7 +65,7 @@ public: bool isMyDomain(int domainIndex) const; // Return processor id where the domain with domainIndex resides - int getProccessorID(int domainIndex) const; + int getProcessorID(int domainIndex) const; //!< Set nb of required domains. (Used to sort joints via jointId()) @@ -127,6 +127,17 @@ const std::vector& domain_meshes); // Evaluate current memory usage and return the maximal one in KB int evaluateMemory() const; + void sendMesh(const ParaMEDMEM::MEDCouplingUMesh& mesh, int target) const; + + void recvMesh(ParaMEDMEM::MEDCouplingUMesh*& mesh, int source) const; + + void sendDoubleVec(const std::vector& vec, int target) const; + + void recvDoubleVec(std::vector& vec, int source) const; + + void sendIntVec(const std::vector& vec, int target) const; + void recvIntVec(std::vector& vec, int source) const; + private: int _rank, _world_size; // my rank and nb of processors diff --git a/src/MEDPartitioner/Makefile.am b/src/MEDPartitioner/Makefile.am index 0032fc256..44bb7573e 100644 --- a/src/MEDPartitioner/Makefile.am +++ b/src/MEDPartitioner/Makefile.am @@ -39,6 +39,7 @@ MEDPARTITIONER_MESHCollectionMedXMLDriver.hxx \ MEDPARTITIONER_MESHCollectionMedAsciiDriver.hxx \ MEDPARTITIONER_ParallelTopology.hxx \ MEDPARTITIONER_FaceModel.hxx \ +MEDPARTITIONER_JointFinder.hxx \ MEDPARTITIONER_Graph.hxx\ MEDPARTITIONER_UserGraph.hxx\ MEDPARTITIONER_SequentialTopology.hxx \ @@ -65,6 +66,7 @@ MEDPARTITIONER_Graph.cxx\ MEDPARTITIONER_UserGraph.cxx\ MEDPARTITIONER_ParaDomainSelector.cxx \ MEDPARTITIONER_JointExchangeData.cxx \ +MEDPARTITIONER_JointFinder.cxx \ MEDPARTITIONER_SkyLineArray.cxx \ MEDPARTITIONER_ConnectZone.cxx diff --git a/src/ParaMEDMEM/ElementLocator.cxx b/src/ParaMEDMEM/ElementLocator.cxx index 3252279cd..e14cf117c 100644 --- a/src/ParaMEDMEM/ElementLocator.cxx +++ b/src/ParaMEDMEM/ElementLocator.cxx @@ -258,7 +258,7 @@ namespace ParaMEDMEM //serialization of local mesh to send data to distant proc. local_mesh->serialize(v1Local,v2Local); //Building the right instance of copy of distant mesh. - MEDCouplingPointSet *distant_mesh_tmp=MEDCouplingPointSet::buildInstanceFromMeshType((MEDCouplingMeshType)tinyInfoDistant[0]); + MEDCouplingPointSet *distant_mesh_tmp=MEDCouplingPointSet::BuildInstanceFromMeshType((MEDCouplingMeshType)tinyInfoDistant[0]); std::vector unusedTinyDistantSts; distant_mesh_tmp->resizeForUnserialization(tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts); int nbLocalElems=0; diff --git a/src/ParaMEDMEM/ICoCoMEDField.cxx b/src/ParaMEDMEM/ICoCoMEDField.cxx index cc3b61dd5..f3f9e66a0 100644 --- a/src/ParaMEDMEM/ICoCoMEDField.cxx +++ b/src/ParaMEDMEM/ICoCoMEDField.cxx @@ -88,6 +88,7 @@ namespace ICoCo //creating a connectivity table that complies to MED (1 indexing) //and passing it to _mesh int* conn=new int[triofield._nodes_per_elem]; + _mesh->setMeshDimension(triofield._mesh_dim); for (int i=0; isetMeshDimension(triofield._mesh_dim); _mesh->finishInsertingCells(); //field on the sending end diff --git a/src/RENUMBER/Makefile.am b/src/RENUMBER/Makefile.am index 397ed5c44..551a5d1f9 100644 --- a/src/RENUMBER/Makefile.am +++ b/src/RENUMBER/Makefile.am @@ -53,7 +53,8 @@ endif librenumber_la_CPPFLAGS= $(MED2_INCLUDES) $(HDF5_INCLUDES) @CXXTMPDPTHFLAGS@ \ $(BOOST_CPPFLAGS) \ - -I$(srcdir)/../MEDMEM -I$(srcdir)/../MEDWrapper/V2_1/Core -I$(srcdir)/../INTERP_KERNEL/Bases + -I$(srcdir)/../MEDMEM -I$(srcdir)/../MEDWrapper/V2_1/Core -I$(srcdir)/../INTERP_KERNEL/Bases \ + -I$(srcdir)/../INTERP_KERNEL/GaussPoints librenumber_la_LDFLAGS= #libmedsplitter_la_LDFLAGS= $(MED2_LIBS) $(HDF5_LIBS) $(STDLIB) \