typedef enum
{
- NORM_POINT0 = 0,
+ NORM_POINT1 = 0,
NORM_SEG2 = 1,
NORM_SEG3 = 2,
NORM_TRI3 = 3,
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
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)));
_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:
#include <sstream>
#include <algorithm>
+#ifdef _POSIX_MAPPED_FILES
+#include <sys/mman.h>
+#else
+#ifdef WNT
+#include <windows.h>
+#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<char> INTERP_KERNEL::AsmX86::convertIntoMachineLangage(const std::vector<std::string>& asmb) const throw(INTERP_KERNEL::Exception)
return ret;
}
-char *INTERP_KERNEL::AsmX86::convertMachineLangageInBasic(const std::vector<char>& ml, int& lgth) const
+char *INTERP_KERNEL::AsmX86::copyToExecMemZone(const std::vector<char>& 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;
}
{
public:
std::vector<char> convertIntoMachineLangage(const std::vector<std::string>& asmb) const throw(INTERP_KERNEL::Exception);
- char *convertMachineLangageInBasic(const std::vector<char>& ml, int& lgth) const;
+ char *copyToExecMemZone(const std::vector<char>& ml, unsigned& offset) const;
private:
void convertOneInstructionInML(const std::string& inst, std::vector<char>& ml) const throw(INTERP_KERNEL::Exception);
private:
#include <iostream>
#include <algorithm>
-#ifdef _POSIX_MAPPED_FILES
-#include <sys/mman.h>
-#else
-#ifdef WNT
-#include <windows.h>
-#endif
-#endif
-
using namespace INTERP_KERNEL;
const char LeafExprVar::END_OF_RECOGNIZED_VAR[]="Vec";
for(std::vector<char>::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
for(std::vector<char>::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<std::string>& ass) const
{
for(std::list<ExprParser>::const_iterator iter=_sub_expr.begin();iter!=_sub_expr.end();iter++)
(*iter).compileX86LowLev(ass);
- for(std::list<Function *>::const_iterator iter2=_func_btw_sub_expr.begin();iter2!=_func_btw_sub_expr.end();iter2++)
- (*iter2)->operateX86(ass);
}
+ for(std::list<Function *>::const_iterator iter2=_func_btw_sub_expr.begin();iter2!=_func_btw_sub_expr.end();iter2++)
+ (*iter2)->operateX86(ass);
}
void ExprParser::compileX86_64LowLev(std::vector<std::string>& ass) const
{
for(std::list<ExprParser>::const_iterator iter=_sub_expr.begin();iter!=_sub_expr.end();iter++)
(*iter).compileX86_64LowLev(ass);
- for(std::list<Function *>::const_iterator iter2=_func_btw_sub_expr.begin();iter2!=_func_btw_sub_expr.end();iter2++)
- (*iter2)->operateX86(ass);
}
+ for(std::list<Function *>::const_iterator iter2=_func_btw_sub_expr.begin();iter2!=_func_btw_sub_expr.end();iter2++)
+ (*iter2)->operateX86(ass);
}
void LeafExprVal::compileX86(std::vector<std::string>& ass) const
--- /dev/null
+// 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 <math.h>
+#include <algorithm>
+#include <sstream>
+
+//#define MYDEBUG
+
+#ifdef MYDEBUG
+#include <iostream>
+#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<<"#######################################################"<<std::endl;
+ std::cout<<"Cell type : "<<info->GetCellType()<<std::endl;
+#endif
+
+ for( int gaussId = 0; gaussId < info->GetNbGauss() ; gaussId++ ) {
+#ifdef MYDEBUG
+ std::cout<<"Gauss ID = "<<gaussId<<std::endl;
+#endif
+
+ double* coord = &aCoords[ gaussId*theSpaceDim ];
+ const double* function = info->GetFunctionValues(gaussId);
+ for ( int connId = 0; connId < aConn ; connId++ ) {
+ const double* nodeCoord = theNodeCoords + (theIndex[connId]*theSpaceDim);
+
+#ifdef MYDEBUG
+ std::cout<<"Function Value["<<connId<<"] = "<<function[connId]<<std::endl;
+
+ std::cout<<"Node #"<<connId <<" ";
+ for( int dimId = 0; dimId < theSpaceDim; dimId++ ) {
+ switch(dimId) {
+ case 0: std::cout<<"( "<<nodeCoord[dimId];break;
+ case 1: std::cout<<", "<<nodeCoord[dimId];break;
+ case 2: std::cout<<", "<<nodeCoord[dimId]<<" )";break;
+ }
+ }
+ std::cout<<std::endl;
+#endif
+
+ for( int dimId = 0; dimId < theSpaceDim; dimId++ ) {
+ coord[dimId] += nodeCoord[dimId]*function[connId];
+ }
+ }
+ }
+ return aCoords;
+}
--- /dev/null
+// 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
+//
+
+#ifndef __INTERPKERNELGAUSS_HXX__
+#define __INTERPKERNELGAUSS_HXX__
+
+#include "NormalizedUnstructuredMesh.hxx"
+#include "InterpKernelException.hxx"
+#include <vector>
+
+namespace INTERP_KERNEL {
+
+ typedef std::vector<double> DataVector;
+ typedef std::vector<int> 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<GaussInfo*> GaussInfoVector;
+ GaussInfoVector myGaussInfo;
+ };
+}
+#endif //INTERPKERNELGAUSS
--- /dev/null
+# 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
#
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
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)
-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 = \
-// 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
%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;
%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;
%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;
%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;
%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;
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
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<int>& tinyInfo) const throw(INTERP_KERNEL::Exception) = 0;
//! size of returned tinyInfo must be always the same.
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<sz;i++)
PyList_SetItem(coords,i,PyFloat_FromDouble(coo[i]));
delete [] coo;
delete [] c;
}
- static void rotate3DAlg(PyObject *center, PyObject *vect, double angle, int nbNodes, PyObject *coords) throw(INTERP_KERNEL::Exception)
+ static void Rotate3DAlg(PyObject *center, PyObject *vect, double angle, int nbNodes, PyObject *coords) throw(INTERP_KERNEL::Exception)
{
int sz,sz2;
double *c=convertPyToNewDblArr2(center,&sz);
double *coo=convertPyToNewDblArr2(coords,&sz);
double *v=convertPyToNewDblArr2(vect,&sz2);
- ParaMEDMEM::MEDCouplingPointSet::rotate3DAlg(c,v,angle,nbNodes,coo);
+ ParaMEDMEM::MEDCouplingPointSet::Rotate3DAlg(c,v,angle,nbNodes,coo);
for(int i=0;i<sz;i++)
PyList_SetItem(coords,i,PyFloat_FromDouble(coo[i]));
delete [] coo;
MEDCouplingFieldDouble *getAspectRatioField() const throw(INTERP_KERNEL::Exception);
MEDCouplingFieldDouble *getWarpField() const throw(INTERP_KERNEL::Exception);
MEDCouplingFieldDouble *getSkewField() const throw(INTERP_KERNEL::Exception);
- static MEDCouplingUMesh *mergeUMeshes(const MEDCouplingUMesh *mesh1, const MEDCouplingUMesh *mesh2) throw(INTERP_KERNEL::Exception);
- static MEDCouplingUMesh *mergeUMeshesOnSameCoords(const MEDCouplingUMesh *mesh1, const MEDCouplingUMesh *mesh2) throw(INTERP_KERNEL::Exception);
+ static MEDCouplingUMesh *MergeUMeshes(const MEDCouplingUMesh *mesh1, const MEDCouplingUMesh *mesh2) throw(INTERP_KERNEL::Exception);
+ static MEDCouplingUMesh *MergeUMeshesOnSameCoords(const MEDCouplingUMesh *mesh1, const MEDCouplingUMesh *mesh2) throw(INTERP_KERNEL::Exception);
%extend {
std::string __str__() const
{
return convertIntArrToPyList2(elts);
}
- static PyObject *mergeUMeshesOnSameCoords(PyObject *ms) throw(INTERP_KERNEL::Exception)
+ static PyObject *MergeUMeshesOnSameCoords(PyObject *ms) throw(INTERP_KERNEL::Exception)
{
- std::vector<ParaMEDMEM::MEDCouplingUMesh *> meshes;
- if(PyList_Check(ms))
- {
- int sz=PyList_Size(ms);
- meshes.resize(sz);
- for(int i=0;i<sz;i++)
- {
- PyObject *o=PyList_GetItem(ms,i);
- void *arg;
- SWIG_ConvertPtr(o,&arg,SWIGTYPE_p_ParaMEDMEM__MEDCouplingUMesh, 0 | 0 );
- meshes[i]=reinterpret_cast<ParaMEDMEM::MEDCouplingUMesh *>(arg);
- }
- }
- else
- {
- PyErr_SetString(PyExc_TypeError,"mergeUMeshesOnSameCoords : not a list as first parameter");
- PyErr_Print();
- return 0;
- }
- MEDCouplingUMesh *ret=MEDCouplingUMesh::mergeUMeshesOnSameCoords(meshes);
+ std::vector<const ParaMEDMEM::MEDCouplingUMesh *> 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<MEDCouplingUMesh *> meshes;
- convertPyObjToVecUMeshes(ms,meshes);
+ std::vector<const MEDCouplingUMesh *> meshes;
+ convertPyObjToVecUMeshesCst(ms,meshes);
std::vector<DataArrayInt *> 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<sz;i++)
return convertDblArrToPyListOfTuple(vals,3,2);
}
- static MEDCouplingUMesh *mergeUMeshes(PyObject *li) throw(INTERP_KERNEL::Exception)
+ static MEDCouplingUMesh *MergeUMeshes(PyObject *li) throw(INTERP_KERNEL::Exception)
{
std::vector<const ParaMEDMEM::MEDCouplingUMesh *> tmp;
convertPyObjToVecUMeshesCst(li,tmp);
- return MEDCouplingUMesh::mergeUMeshes(tmp);
+ return MEDCouplingUMesh::MergeUMeshes(tmp);
}
PyObject *areCellsIncludedIn(const MEDCouplingUMesh *other, int compType) const throw(INTERP_KERNEL::Exception)
}
void convertToPolyTypes(const std::vector<int>& 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
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);
}
return convertDblArrToPyList(tmp,sz);
}
- static DataArrayDouble *aggregate(PyObject *li) throw(INTERP_KERNEL::Exception)
+ static DataArrayDouble *Aggregate(PyObject *li) throw(INTERP_KERNEL::Exception)
{
std::vector<const DataArrayDouble *> 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<const DataArrayDouble *> tmp;
convertPyObjToVecDataArrayDblCst(li,tmp);
- return DataArrayDouble::meld(tmp);
+ return DataArrayDouble::Meld(tmp);
}
};
return self->repr();
}
+ PyObject *getDifferentValues(bool val) const throw(INTERP_KERNEL::Exception)
+ {
+ std::set<int> 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);
}
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<DataArrayInt *> groups;
+ std::vector<const DataArrayInt *> groups;
std::vector< std::vector<int> > 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();
return convertIntArrToPyList(tmp,sz);
}
- static DataArrayInt *meld(PyObject *li) throw(INTERP_KERNEL::Exception)
+ static DataArrayInt *Meld(PyObject *li) throw(INTERP_KERNEL::Exception)
+ {
+ std::vector<const DataArrayInt *> tmp;
+ convertPyObjToVecDataArrayIntCst(li,tmp);
+ return DataArrayInt::Meld(tmp);
+ }
+
+ static DataArrayInt *BuildUnion(PyObject *li) throw(INTERP_KERNEL::Exception)
+ {
+ std::vector<const DataArrayInt *> tmp;
+ convertPyObjToVecDataArrayIntCst(li,tmp);
+ return DataArrayInt::BuildUnion(tmp);
+ }
+
+ static DataArrayInt *BuildIntersection(PyObject *li) throw(INTERP_KERNEL::Exception)
{
std::vector<const DataArrayInt *> tmp;
convertPyObjToVecDataArrayIntCst(li,tmp);
- return DataArrayInt::meld(tmp);
+ return DataArrayInt::BuildIntersection(tmp);
}
PyObject *getMaxValue() const throw(INTERP_KERNEL::Exception)
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);
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);
self->setSelectedComponents(f,tmp);
}
- static MEDCouplingFieldDouble *mergeFields(PyObject *li) throw(INTERP_KERNEL::Exception)
+ static MEDCouplingFieldDouble *MergeFields(PyObject *li) throw(INTERP_KERNEL::Exception)
{
std::vector<const MEDCouplingFieldDouble *> tmp;
convertPyObjToVecFieldDblCst(li,tmp);
- return MEDCouplingFieldDouble::mergeFields(tmp);
+ return MEDCouplingFieldDouble::MergeFields(tmp);
}
}
};
#endif
}
+static PyObject *convertIntArrToPyList3(const std::set<int>& v) throw(INTERP_KERNEL::Exception)
+{
+ int size=v.size();
+ PyObject *ret=PyList_New(size);
+ std::set<int>::const_iterator it=v.begin();
+ for(int i=0;i<size;i++,it++)
+ PyList_SetItem(ret,i,PyInt_FromLong(*it));
+ return ret;
+}
+
static PyObject *convertIntArrToPyListOfTuple(const int *vals, int nbOfComp, int nbOfTuples) throw(INTERP_KERNEL::Exception)
{
PyObject *ret=PyList_New(nbOfTuples);
}
}
+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<size;i++)
+ {
+ PyObject *o=PyList_GetItem(pyLi,i);
+ if(PyInt_Check(o))
+ {
+ int val=(int)PyInt_AS_LONG(o);
+ if(i<sizeOfArray)
+ arrToFill[i]=val;
+ }
+ else
+ {
+ const char msg[]="list must contain integers only";
+ PyErr_SetString(PyExc_TypeError,msg);
+ throw INTERP_KERNEL::Exception(msg);
+ }
+ }
+ for(int i=size;i<sizeOfArray;i++)
+ arrToFill[i]=dftVal;
+ return;
+
+ }
+ else if(PyTuple_Check(pyLi))
+ {
+ int size=PyTuple_Size(pyLi);
+ for(int i=0;i<size;i++)
+ {
+ PyObject *o=PyTuple_GetItem(pyLi,i);
+ if(PyInt_Check(o))
+ {
+ int val=(int)PyInt_AS_LONG(o);
+ if(i<sizeOfArray)
+ arrToFill[i]=val;
+ }
+ else
+ {
+ const char msg[]="tuple must contain integers only";
+ PyErr_SetString(PyExc_TypeError,msg);
+ throw INTERP_KERNEL::Exception(msg);
+ }
+ }
+ for(int i=size;i<sizeOfArray;i++)
+ arrToFill[i]=dftVal;
+ return;
+ }
+ else
+ {
+ const char msg[]="fillArrayWithPyListInt : not a list";
+ PyErr_SetString(PyExc_TypeError,msg);
+ throw INTERP_KERNEL::Exception(msg);
+ }
+}
+
static PyObject *convertDblArrToPyList(const double *ptr, int size) throw(INTERP_KERNEL::Exception)
{
PyObject *ret=PyList_New(size);
}
else
{
- const char msg[]="convertPyToNewIntArr : not a list";
+ const char msg[]="convertPyToNewDblArr2 : not a list";
PyErr_SetString(PyExc_TypeError,msg);
throw INTERP_KERNEL::Exception(msg);
}
}
-void convertPyObjToVecUMeshes(PyObject *ms, std::vector<ParaMEDMEM::MEDCouplingUMesh *>& 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<size;i++)
{
- PyObject *obj=PyList_GetItem(ms,i);
- void *argp;
- int status=SWIG_ConvertPtr(obj,&argp,SWIGTYPE_p_ParaMEDMEM__MEDCouplingUMesh,0|0);
- if(!SWIG_IsOK(status))
+ PyObject *o=PyList_GetItem(pyLi,i);
+ if(PyFloat_Check(o))
{
- const char msg[]="list must contain only MEDCouplingUMesh";
+ double val=PyFloat_AS_DOUBLE(o);
+ if(i<sizeOfArray)
+ arrToFill[i]=val;
+ }
+ else if(PyInt_Check(o))
+ {
+ long val0=PyInt_AS_LONG(o);
+ double val=val0;
+ if(i<sizeOfArray)
+ arrToFill[i]=val;
+ }
+ else
+ {
+ const char msg[]="fillArrayWithPyListDbl : list must contain floats/integers only";
+ PyErr_SetString(PyExc_TypeError,msg);
+ throw INTERP_KERNEL::Exception(msg);
+ }
+ }
+ for(int i=size;i<sizeOfArray;i++)
+ arrToFill[i]=dftVal;
+ return;
+ }
+ else if(PyTuple_Check(pyLi))
+ {
+ int size=PyTuple_Size(pyLi);
+ for(int i=0;i<size;i++)
+ {
+ PyObject *o=PyTuple_GetItem(pyLi,i);
+ if(PyFloat_Check(o))
+ {
+ double val=PyFloat_AS_DOUBLE(o);
+ arrToFill[i]=val;
+ }
+ else if(PyInt_Check(o))
+ {
+ long val0=PyInt_AS_LONG(o);
+ double val=val0;
+ arrToFill[i]=val;
+ }
+ else
+ {
+ const char msg[]="fillArrayWithPyListDbl : tuple must contain floats/integers only";
PyErr_SetString(PyExc_TypeError,msg);
throw INTERP_KERNEL::Exception(msg);
}
- ParaMEDMEM::MEDCouplingUMesh *arg=reinterpret_cast< ParaMEDMEM::MEDCouplingUMesh * >(argp);
- v[i]=arg;
}
+ for(int i=size;i<sizeOfArray;i++)
+ arrToFill[i]=dftVal;
+ return ;
}
else
{
- const char msg[]="convertPyObjToVecUMeshes : not a list";
+ const char msg[]="convertPyToNewIntArr : not a list";
PyErr_SetString(PyExc_TypeError,msg);
throw INTERP_KERNEL::Exception(msg);
}
}
}
-void convertPyObjToVecDataArrayInt(PyObject *ms, std::vector<ParaMEDMEM::DataArrayInt *>& v) throw(INTERP_KERNEL::Exception)
-{
- if(PyList_Check(ms))
- {
- int size=PyList_Size(ms);
- v.resize(size);
- for(int i=0;i<size;i++)
- {
- PyObject *obj=PyList_GetItem(ms,i);
- void *argp;
- int status=SWIG_ConvertPtr(obj,&argp,SWIGTYPE_p_ParaMEDMEM__DataArrayInt,0|0);
- if(!SWIG_IsOK(status))
- {
- const char msg[]="list must contain only DataArrayInt";
- PyErr_SetString(PyExc_TypeError,msg);
- throw INTERP_KERNEL::Exception(msg);
- }
- ParaMEDMEM::DataArrayInt *arg=reinterpret_cast< ParaMEDMEM::DataArrayInt * >(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<const ParaMEDMEM::DataArrayInt *>& v) throw(INTERP_KERNEL::Exception)
{
if(PyList_Check(ms))
}
+void MEDFileUMesh::copyFamGrpMapsFrom(const MEDFileUMesh& other)
+{
+ _groups=other._groups;
+ _families=other._families;
+}
+
+void MEDFileUMesh::setFamilyInfo(const std::map<std::string,int>& info)
+{
+ _families=info;
+}
+
+void MEDFileUMesh::setGroupInfo(const std::map<std::string, std::vector<std::string> >&info)
+{
+ _groups=info;
+}
+
void MEDFileUMesh::setFamilyNameAttachedOnId(int id, const std::string& newFamName) throw(INTERP_KERNEL::Exception)
{
std::string oldName=getFamilyNameGivenId(id);
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<std::string,int>& getFamilyInfo() const { return _families; }
+ const std::map<std::string, std::vector<std::string> >& 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<std::string,int>& info);
+ void setGroupInfo(const std::map<std::string, std::vector<std::string> >&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<const DataArrayInt *>& grps, bool renum=true) throw(INTERP_KERNEL::Exception);
%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
{
}
static void WriteUMeshesPartition(const char *fileName, const char *meshName, PyObject *li, bool writeFromScratch) throw(INTERP_KERNEL::Exception)
{
- std::vector<ParaMEDMEM::MEDCouplingUMesh *> v=convertFieldDoubleVecFromPy(li);
+ std::vector<const ParaMEDMEM::MEDCouplingUMesh *> 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<ParaMEDMEM::MEDCouplingUMesh *> v=convertFieldDoubleVecFromPy(li);
+ std::vector<const ParaMEDMEM::MEDCouplingUMesh *> 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<ParaMEDMEM::MEDCouplingUMesh *> v=convertFieldDoubleVecFromPy(li);
- std::vector<const ParaMEDMEM::MEDCouplingUMesh *> v2(v.begin(),v.end());
- MEDLoader::WriteUMeshes(fileName,v2,writeFromScratch);
+ std::vector<const ParaMEDMEM::MEDCouplingUMesh *> v=convertUMeshVecFromPy(li);
+ MEDLoader::WriteUMeshes(fileName,v,writeFromScratch);
}
static PyObject *GetTypesOfField(const char *fileName, const char *fieldName, const char *meshName) throw(INTERP_KERNEL::Exception)
{
%include "MEDFileMesh.hxx"
+%extend ParaMEDMEM::MEDFileUMesh
+{
+ void setGroupsAtLevel(int meshDimRelToMaxExt, PyObject *li, bool renum=true) throw(INTERP_KERNEL::Exception)
+ {
+ std::vector<const DataArrayInt *> grps;
+ convertPyObjToVecDataArrayIntCst(li,grps);
+ self->setGroupsAtLevel(meshDimRelToMaxExt,grps,renum);
+ }
+}
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
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")
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
#include "MEDCouplingNormalizedUnstructuredMesh.hxx"
#include "MEDCouplingMemArray.hxx"
#include "PointLocator3DIntersectorP0P0.hxx"
-
+#include "BBTree.txx"
#include "MEDPARTITIONER_utils.hxx"
#include "MEDPARTITIONER_Graph.hxx"
#include "MEDPARTITIONER_MESHCollectionDriver.hxx"
#include "MEDPARTITIONER_MESHCollectionMedXMLDriver.hxx"
#include "MEDPARTITIONER_MESHCollectionMedAsciiDriver.hxx"
-
+#include "MEDPARTITIONER_JointFinder.hxx"
#include "MEDPARTITIONER_UserGraph.hxx"
#include <vector>
#include <string>
+#include <limits>
#ifndef WNT
# include <ext/hash_map>
_create_empty_groups(create_empty_groups)
{
- _mesh.resize(_topology->nbDomain());
-
- //splitting the initial domains into smaller bits
-
- std::vector<std::vector<ParaMEDMEM::MEDCouplingUMesh*> > splitMeshes;
- splitMeshes.resize(topology->nbDomain());
- for (int inew=0; inew<topology->nbDomain();inew++)
- splitMeshes[inew].resize(initialCollection.getTopology()->nbDomain());
-
std::vector<std::vector<std::vector<int> > > new2oldIds(initialCollection.getTopology()->nbDomain());
+ castCellMeshes(initialCollection, new2oldIds);
- for (int iold=0; iold<initialCollection.getTopology()->nbDomain();iold++)
- {
- if (!isParallelMode() || initialCollection._domain_selector->isMyDomain(iold))
- {
- int size=(initialCollection._mesh)[iold]->getNumberOfCells();
- std::vector<int> globalids(size);
- initialCollection.getTopology()->getCellList(iold, &globalids[0]);
- std::vector<int> ilocalnew(size);
- std::vector<int> ipnew(size);
- topology->convertGlobalCellList(&globalids[0],size,&ilocalnew[0],&ipnew[0]);
- new2oldIds[iold].resize(topology->nbDomain());
- for (int i=0; i<ilocalnew.size();i++)
- {
- new2oldIds[iold][ipnew[i]].push_back(i);
- }
- for (int inew=0;inew<topology->nbDomain();inew++)
- {
- // cout <<"inew:"<<inew<<" iold:"<<iold<<endl;
- // for (int i=0; i<new2oldIds[iold][inew].size();i++)
- // cout<<new2oldIds[iold][inew][i]<<" ";
- // cout<<endl;
- splitMeshes[inew][iold]=(ParaMEDMEM::MEDCouplingUMesh*)(initialCollection.getMesh())[iold]->buildPartOfMySelf(&new2oldIds[iold][inew][0],&new2oldIds[iold][inew][0]+new2oldIds[iold][inew].size(),true);
- // cout <<"small domain "<<iold<<" "<<inew<<" has "<<splitMeshes[inew][iold]->getNumberOfCells()<<" cells"<<endl;
- }
- }
- }
- if (isParallelMode())
- {
- //send/receive stuff
- }
-
- //fusing the split meshes
- for (int inew=0; inew<topology->nbDomain() ;inew++)
- {
- std::vector<const ParaMEDMEM::MEDCouplingUMesh*> 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 #"<<inew+1<<" with "<<_mesh[inew]->getNumberOfCells()<<" cells and "<<_mesh[inew]->getNumberOfNodes()<<" nodes"<<std::endl;
- }
- for (int i=0; i< splitMeshes[inew].size();i++)
- splitMeshes[inew][i]->decrRef();
- }
//casting cell families on new meshes
_cellFamilyIds.resize(topology->nbDomain());
castIntField(initialCollection.getMesh(), this->getMesh(),initialCollection.getCellFamilyIds(),_cellFamilyIds, new2oldIds);
// treating faces
/////////////////
- std::multimap<pair<int,int>, pair<int,int> > nodeMapping;
+ NodeMapping nodeMapping;
createNodeMapping(initialCollection, nodeMapping);
std::vector<std::vector<std::vector<int> > > new2oldFaceIds;
castMeshes(initialCollection.getFaceMesh(), this->getFaceMesh(),initialCollection, nodeMapping, new2oldFaceIds);
+
+ ////////////////////
+ //treating families
+ ////////////////////
+
_faceFamilyIds.resize(topology->nbDomain());
//allocating family ids arrays
for (int inew=0; inew<topology->nbDomain();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);
}
}
+/*!
+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<std::vector<std::vector<int> > >& 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<std::vector<ParaMEDMEM::MEDCouplingUMesh*> > 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; iold<initialCollection.getTopology()->nbDomain();iold++)
+ {
+ if (!isParallelMode() || initialCollection._domain_selector->isMyDomain(iold))
+ {
+ int size=(initialCollection._mesh)[iold]->getNumberOfCells();
+ std::vector<int> globalids(size);
+ initialCollection.getTopology()->getCellList(iold, &globalids[0]);
+ std::vector<int> ilocalnew(size);
+ std::vector<int> ipnew(size);
+ _topology->convertGlobalCellList(&globalids[0],size,&ilocalnew[0],&ipnew[0]);
+ new2oldIds[iold].resize(_topology->nbDomain());
+ for (int i=0; i<ilocalnew.size();i++)
+ {
+ new2oldIds[iold][ipnew[i]].push_back(i);
+ }
+ for (int inew=0;inew<_topology->nbDomain();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; iold<initialCollection.getTopology()->nbDomain();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"<<iold<<" inew "<<inew<<" "<<splitMeshes[inew][iold]->getNumberOfCells()<<std::endl;
+
+ }
+ if (_domain_selector->isMyDomain(inew))
+ {
+ _domain_selector->recvMesh(splitMeshes[inew][iold],_domain_selector->getProcessorID(iold));
+ std::cout<<"recv iold"<<iold<<" inew "<<inew<<std::endl;
+ }
+ }
+ }
+
+ //fusing the split meshes
+ for (int inew=0; inew<_topology->nbDomain() ;inew++)
+ {
+ std::vector<const ParaMEDMEM::MEDCouplingUMesh*> 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"<<meshes.size()<<std::endl;
+
+ if (!isParallelMode()||_domain_selector->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<int,int>,pair<int,int> >& nodeMapping)
+void MESHCollection::createNodeMapping( MESHCollection& initialCollection, NodeMapping& nodeMapping)
{
+
+ // NodeMapping reverseNodeMapping;
for (int iold=0; iold<initialCollection.getTopology()->nbDomain();iold++)
{
- std::map<pair<double,pair<double, double> >, int > nodeClassifier;
- for (int inode=0; inode<initialCollection.getMesh(iold)->getNumberOfNodes(); inode++)
+
+ double* bbox;
+ BBTree<3>* tree;
+ if (!isParallelMode() || (_domain_selector->isMyDomain(iold)))
{
+ // std::map<pair<double,pair<double, double> >, 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<double, pair<double,double> > 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<nvertices*3;i++)
+ {
+ bbox[i*2]=coordsPtr[i]-1e-9;
+ bbox[i*2+1]=coordsPtr[i]+1e-9;
+ }
+ tree=new BBTree<3>(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<pair<double,pair<double, double> >, 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 <<"("<<iold<<","<<iter->second<<")-->("<<inew<<","<<inode<<")"<<endl;
+ std::cout<<"sendTo"<<_domain_selector->getProcessorID(iold)<<std::endl;
+ _domain_selector->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) <<std::endl;
+ _domain_selector->recvMesh(mesh, _domain_selector->getProcessorID(inew));
+ ParaMEDMEM::DataArrayDouble* coords = mesh->getCoords();
+ for (int inode=0; inode<mesh->getNumberOfNodes();inode++)
+ {
+ double* coordsPtr=coords->getPointer()+inode*3;
+ std::vector<int> 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<int> 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 :" <<inode<<" ("<<iold<<","<<elems[0]<<")-->("<<inew<<","<<inode<<")"<<endl;
+ }
}
}
}
+ std::cout<<"NodeMapping size"<<nodeMapping.size()<<std::endl;
+
}
/*!
faces at the interface are duplicated
*/
-void MESHCollection::castMeshes(std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastFrom,std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastTo, MESHCollection& initialCollection,const multimap<pair<int,int>,pair<int,int> >& nodeMapping, std::vector<std::vector<std::vector<int> > >& new2oldIds)
+void MESHCollection::castMeshes(std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastFrom,std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastTo, MESHCollection& initialCollection,const NodeMapping& nodeMapping, std::vector<std::vector<std::vector<int> > >& new2oldIds)
{
//splitMeshes structure will contain the partition of
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; iold<meshesCastFrom.size();iold++)
{
+ if (isParallelMode() && !_domain_selector->isMyDomain(iold)) continue;
new2oldIds[iold].resize(newSize);
for (int ielem=0;ielem<meshesCastFrom[iold]->getNumberOfCells();ielem++)
{
for (int inode=0;inode<nodes.size();inode++)
{
typedef multimap<pair<int,int>,pair<int,int> >::const_iterator MI;
- pair <MI,MI> myRange = nodeMapping.equal_range(make_pair(iold,nodes[inode]));
+ int mynode=nodes[inode];
+ if (mynode <0 || mynode > 1000000000) exit(1);
+ pair <MI,MI> myRange = nodeMapping.equal_range(make_pair(iold,mynode));
// cout << iold <<" " <<nodes[inode]<<endl;
for (MI iter=myRange.first; iter!=myRange.second; iter++)
{
//creating the splitMeshes from the face ids
for (int inew=0;inew<_topology->nbDomain();inew++)
{
- cout<<"nb faces "<<new2oldIds[iold][inew].size()<<endl;
+ cout<<"nb faces - iold "<<iold<<" inew "<<inew<<" : "<<new2oldIds[iold][inew].size()<<endl;
splitMeshes[inew][iold]=(ParaMEDMEM::MEDCouplingUMesh*)(meshesCastFrom[iold]->buildPartOfMySelf(&new2oldIds[iold][inew][0],&new2oldIds[iold][inew][0]+new2oldIds[iold][inew].size(),true));
}
}
// send/receive stuff
+ if (isParallelMode())
+ for (int iold=0; iold<meshesCastFrom.size();iold++)
+ for (int inew=0; inew<newSize; inew++)
+ {
+ if (_domain_selector->isMyDomain(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<const ParaMEDMEM::MEDCouplingUMesh*> myMeshes(meshesCastFrom.size());
+ vector<const ParaMEDMEM::MEDCouplingUMesh*> myMeshes;
for (int iold=0; iold < meshesCastFrom.size();iold++)
{
- myMeshes[iold]=splitMeshes[inew][iold];
- cout<<"number of nodes"<<splitMeshes[inew][iold]->getNumberOfNodes()<<endl;
- cout<<"number of cells"<<splitMeshes[inew][iold]->getNumberOfCells()<<endl;
- cout<<"dimension "<<splitMeshes[inew][iold]->getMeshDimension()<<endl;
-
+ if (splitMeshes[inew][iold] !=0)
+ myMeshes.push_back(splitMeshes[inew][iold]);
}
meshesCastTo[inew]=ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes(myMeshes);
- meshesCastTo[inew]->checkCoherency();
meshesCastTo[inew]->zipCoords();
- meshesCastTo[inew]->checkCoherency();
for (int iold=0; iold < meshesCastFrom.size();iold++)
- splitMeshes[inew][iold]->decrRef();
- cout<<"creating face mesh #"<<inew<<" with "<<meshesCastTo[inew]->getNumberOfCells()<<" cells and "<<meshesCastTo[inew]->getNumberOfNodes()<<" nodes"<<endl;
+ if (splitMeshes[inew][iold]!=0) splitMeshes[inew][iold]->decrRef();
}
}
vector<const ParaMEDMEM::DataArrayInt*> 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()));
}
{
return _faceMesh;
}
-ParaMEDMEM::MEDCouplingUMesh* MESHCollection::getMesh(int idomain)
+ParaMEDMEM::MEDCouplingUMesh* MESHCollection::getMesh(int idomain) const
{
return _mesh[idomain];
}
multimap< int, int > node2cell;
multimap< int, int > cell2cell;
+ vector<vector<multimap<int,int> > > 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; idomain<nbdomain;idomain++)
{
+ if (isParallelMode() && !_domain_selector->isMyDomain(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);
index_ptr[_mesh[idomain]->getNumberOfNodes()],
globalRevConn->getPointer());
- int min=10000000,max=-1;
+ int min=std::numeric_limits<int>::max(),max=-1;
for (int i=0; i< index_ptr[_mesh[idomain]->getNumberOfNodes()];i++)
{
if (revConn->getPointer()[i]<min) min=revConn->getPointer()[i];
cout <<"min:"<<min<<" max:"<<max<<endl;
int* globalNodeIds=new int[_mesh[idomain]->getNumberOfNodes()];
_topology->getNodeList(idomain,globalNodeIds);
-
+
int* globalRevConnPtr=globalRevConn->getPointer();
for (int i=0; i<_mesh[idomain]->getNumberOfNodes();i++)
{
for (int icell=index_ptr[i]; icell<index_ptr[i+1];icell++)
- node2cell.insert(make_pair(globalNodeIds[i],globalRevConnPtr[icell]));
+ node2cell.insert(make_pair(globalNodeIds[i],globalRevConnPtr[icell]+procOffset));
+ }
+ if (isParallelMode())
+ {
+ for (int mydomain=0; mydomain<_topology->nbDomain();mydomain++)
+ {
+ if (_domain_selector->isMyDomain(mydomain)) continue;
+ // for (int idomain=0; idomain<_topology->nbDomain();idomain++)
+ // {
+ // if (_domain_selector->isMyDomain(idomain)) continue;
+ multimap<int,int>::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 "<<globalNodeIds[ilocnode]<<" "<<icell+offset<<" "<<offset<<endl;
+ }
+ // }
+ }
}
+
globalRevConn->decrRef();
revConn->decrRef();
index->decrRef();
- delete[]globalNodeIds;
+ delete[] globalNodeIds;
}
//creating graph arcs (cell to cell relations)
// 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; i<nbdomain;i++)
+ if (_domain_selector->isMyDomain(i)) maxcell+=_mesh[i]->getNumberOfCells();
+
+ }
+ else
+ {
+ mincell=0;
+ maxcell=_topology->nbCells();
+ }
+ cout<<"mincell"<<mincell<<" maxcell "<<maxcell<<endl;
for (int inode=0; inode<_topology->nbNodes();inode++)
{
typedef multimap<int,int>::const_iterator MI;
{
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->second<maxcell) cell2cell.insert(make_pair(cell1->second,cell2->second));
}
}
}
}
array=new MEDPARTITIONER::MEDSKYLINEARRAY(index,value);
- cout<<"taille du graphe créé "<<array->getNumberOf()<<endl;
- // int cell_number=1;
- // int node_number=1;
- // for (int i=0; i<_topology->nbDomain(); i++)
- // {
- // cell_number+=_topology->getCellNumber(i);
- // node_number+=_topology->getNodeNumber(i);
- // }
- // //list of cells for a given node
- // //vector< vector<int> > node2cell(node_number);
- // map< int, vector<int> > node2cell;
-
- // //list of nodes for a given cell
- // //vector< vector <int> > cell2node(cell_number);
- // map< int, vector <int> > cell2node;
-
- // // map<MED_EN::medGeometryElement,int*> 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"<<endl;
- // //creating the MEDMEMSKYLINEARRAY containing the graph
-
- // int* size = new int[_topology->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<int> 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<int> 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<int,int> cells_neighbours;
- // for (int i=0; i< _topology->nbCells(); i++)
- // {
-
-
- // vector<int> cells(50);
-
- // // /*// cout << "size cell2node "<<cell2node[i+1].size()<<endl;
- // // for (vector<int>::const_iterator iternode=cell2node[i+1].begin();
- // // iternode!=cell2node[i+1].end();
- // // iternode++)
- // // {
- // // int nodeid=*iternode;
- // // // cout << "size node2cell "<<node2cell[nodeid].size()<<endl;
- // // for (vector<int>::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 "<<node2cell[nodeid].size()<<endl;
- // for (int icell=0; icell<node2cell[nodeid].size();icell++)
- // cells_neighbours[node2cell[nodeid][icell]]++;
- // }
- // size[i]=0;
- // int dimension = getMeshDimension();
- // cells.clear();
-
- // for (map<int,int>::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()<<endl;
- // //cout << cells_neighbours.size()<<endl;
-
- // temp[i]=new int[size[i]];
- // // if (has_indivisible_regions)
- // // temp_edgeweight[i]=new int[size[i]];
- // //
- // int itemp=0;
-
- // // memcpy(temp[i],cells,isize*sizeof(int));
-
- // for (vector<int>::const_iterator iter=cells.begin(); iter!=cells.end();iter++)
- // //for(int j=0; j<isize; j++)
- // {
- // temp[i][itemp]=*iter;
- // //temp[i][itemp]=cells[j];
- // // if (has_indivisible_regions)
- // // {
- // // int tag1 = indivisible_tag[(i+1)-1];
- // // //int tag2 = indivisible_tag[iter->first-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"<<endl;
- // int* index=new int[_topology->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<index[i+1];j++)
- // // edgeweights[j-1]=temp_edgeweight[i][j-index[i]];
- // // delete[] temp_edgeweight[i];
- // // }
- // // delete[]temp_edgeweight;
- // // }
- // delete[] index;
- // delete[] temp;
- // delete[] size;
cout<< "end of graph creation"<<endl;
}
+// /*! Method contributing to the distant cell graph
+// */
+// void MESHCollection::findCommonDistantNodes(vector<vector<multimap<int,int> > >& commonDistantNodes)
+// {
+// int nbdomain=_topology->nbDomain();
+// commonDistantNodes.resize(nbdomain);
+// for (int i=0; i<nbdomain;i++) commonDistantNodes[i].resize(nbdomain);
+// int nbproc=_domain_selector->nbProcs();
+// vector<BBTree<3>* > bbtree(nbdomain);
+// vector<ParaMEDMEM::DataArrayInt*> rev(nbdomain);
+// vector<ParaMEDMEM::DataArrayInt*>revIndx(nbdomain);
+// int meshDim;
+// int spaceDim;
+
+// for (int mydomain=0;mydomain<nbdomain;mydomain++)
+// {
+// if(! _domain_selector->isMyDomain(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;isource<nbdomain;isource++)
+// for (int itarget=0;itarget<nbdomain;itarget++)
+// {
+
+// if (_domain_selector->isMyDomain(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<double> 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<int> localCorrespondency;
+// _domain_selector->recvIntVec(localCorrespondency, targetProc);
+// cout<<"size "<<localCorrespondency.size()<<endl;
+// for (int i=0; i<localCorrespondency.size()/2;i++)
+// commonDistantNodes[isource][itarget].insert(make_pair(localCorrespondency[2*i],localCorrespondency[2*i+1]));
+
+// }
+// if (_domain_selector->isMyDomain(itarget))
+// {
+// //receiving data from source proc
+// int sourceProc = isource%nbproc;
+// std::vector<double> recvVec;
+// _domain_selector->recvDoubleVec(recvVec,sourceProc);
+// std::map<int,int> 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<spaceDim;i++)
+// {
+// bbox[2*i]=recvVec[inode*spaceDim+i]-1e-12;
+// bbox[2*i+1]=bbox[2*i]+2e-12;
+// }
+// vector<int> inodes;
+// bbtree[itarget]->getIntersectingElems(bbox,inodes);
+// delete[] bbox;
+
+// if (inodes.size()>0) commonNodes.insert(make_pair(inodes[0],inode));
+// }
+// std::vector<int> nodeCellCorrespondency;
+// for (map<int,int>::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];icell<revIndxPtr[iter->first+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
// }
// }
-// 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<int>(old_collection, field_char, itnumber, ordernumber);
-// else
-// castFields<double>(old_collection, field_char, itnumber, ordernumber);
+// std::vector<ParaMEDMEM::TypeOfField> fieldTypes =MEDLoader::GetTypesOfField(fileName, fieldName,meshName);
+
// }
// void MESHCollection::castAllFields(const MESHCollection& initial_collection)
// vector <int> iternumber;
// vector <int> ordernumber;
// vector <int> 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<field_names.size(); i++)
// {
// char field_char[80];
// strcpy(field_char,field_names[i].c_str());
-// // choosing whether the field is of int or double type
-// if (types[i] ==0)
-// castFields<int>(initial_collection, field_char, iternumber[i], ordernumber[i]);
-// else
-// castFields<double>(initial_collection, field_char, iternumber[i], ordernumber[i]);
+// castFieldDouble(initial_collection, field_char, iternumber[i], ordernumber[i]);
// }
// }
{
ostringstream oss;
oss<<name<<"_"<<i;
+ if (!isParallelMode() || _domain_selector->isMyDomain(i))
_mesh[i]->setName(oss.str().c_str());
}
}
//#include "boost/shared_ptr.hpp"
#include <vector>
#include <map>
-
+#include <string>
namespace ParaMEDMEM
{
class MEDCouplingUMesh;
class ParaDomainSelector;
class MEDSKYLINEARRAY;
class CONNECTZONE;
-
+ class JointFinder;
typedef enum{MedAscii, MedXML, Undefined} DriverType;
+ typedef std::multimap<std::pair<int,int>, std::pair<int,int> > NodeMapping ;
+ typedef std::vector<std::pair<int,int> > NodeList;
class MEDPARTITIONER_EXPORT MESHCollection
{
std::vector<ParaMEDMEM::MEDCouplingUMesh*>& getFaceMesh() ;
std::vector<std::vector<ParaMEDMEM::MEDCouplingUMesh*> >&getGroupMeshes();
- ParaMEDMEM::MEDCouplingUMesh* getMesh(int idomain);
+ ParaMEDMEM::MEDCouplingUMesh* getMesh(int idomain) const;
ParaMEDMEM::MEDCouplingUMesh* getFaceMesh(int idomain);
std::vector<ParaMEDMEM::MEDCouplingUMesh*>& getGroupMeshes(int idomain);
//getting a pointer to topology
Topology* getTopology() const ;
-
+ ParaDomainSelector* getParaDomainSelector() const{return _domain_selector;}
//settig a new topology
void setTopology(Topology* topology);
//creates the node mapping between an old collection and the present one
void createNodeMapping( MESHCollection& initialCollection, std::multimap<std::pair<int,int>,std::pair<int,int> >& nodeMapping);
+ void castCellMeshes(MESHCollection& initialCollection, std::vector<std::vector<std::vector<int> > >& new2oldIds);
//creates faces on the new collection
void castFaceMeshes( MESHCollection& initialCollection,const std::multimap<std::pair<int,int>,std::pair<int,int> >& nodeMapping);
void castIntField(std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastFrom,std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastTo, std::vector<ParaMEDMEM::DataArrayInt*>& arrayFrom, std::vector<ParaMEDMEM::DataArrayInt*>& arrayTo, std::vector<std::vector< std::vector<int> > > & old2newMapping);
+ void findCommonDistantNodes(std::vector<std::vector<std::multimap<int,int> > >& commonDistantNodes);
// template <class T>
// void castFields(const MESHCollection& old_collection, const string& fieldname, int itnumber, int ordernumber);
/*! flag specifying that groups must be created on all domains,
even if they are empty*/
bool _create_empty_groups;
+
+ JointFinder* _joint_finder;
};
}//of namespace
int edgecut;
int* partition = new int[n];
+ cout << "ParMETIS : n="<<n<<endl;
if (nparts >1)
{
if ( parallelizer )
// distribution of vertices of the graph among the processors
int * vtxdist = parallelizer ? parallelizer->getNbVertOfProcs() : 0;
MPI_Comm comm = MPI_COMM_WORLD;
-
+ cout<<"vtxdist[1]"<<" "<<vtxdist[2]<<endl;
ParMETIS_PartKway( vtxdist, xadj, adjncy, vwgt, adjwgt, &wgtflag,
&base, &nparts, options, &edgecut, partition, &comm );
#else
bool ParaDomainSelector::isMyDomain(int domainIndex) const
{
evaluateMemory();
- return (_rank == getProccessorID( domainIndex ));
+ return (_rank == getProcessorID( domainIndex ));
}
//================================================================================
*/
//================================================================================
-int ParaDomainSelector::getProccessorID(int domainIndex) const
+int ParaDomainSelector::getProcessorID(int domainIndex) const
{
evaluateMemory();
return ( domainIndex % _world_size );
ordered_nbs.push_back(0);
for ( int iproc = 0; iproc < nbProcs(); ++iproc )
for ( int idomain = 0; idomain < nb_domains; ++idomain )
- if ( getProccessorID( idomain ) == iproc )
+ if ( getProcessorID( idomain ) == iproc )
{
domain_order[ idomain ] = ordered_nbs.size() - 1;
ordered_nbs.push_back( ordered_nbs.back() + all_nb_elems[idomain] );
elem_shift_by_domain.resize( nb_domains+1, 0 );
for ( int idomain = 0; idomain < nb_domains; ++idomain )
elem_shift_by_domain[ idomain ] = ordered_nbs[ domain_order[ idomain ]];
-
+
elem_shift_by_domain.back() = ordered_nbs.back(); // to know total nb of elements
-
- // if ( entity == MED_CELL )
- {
- // fill _nb_vert_of_procs
- _nb_vert_of_procs.resize( _world_size+1, 0 );
- for ( int i = 0; i < nb_domains; ++i )
+
+ // fill _nb_vert_of_procs
+ _nb_vert_of_procs.resize( _world_size+1, 0 );
+ for ( int i = 0; i < nb_domains; ++i )
{
- int rank = getProccessorID( i );
+ int rank = getProcessorID( i );
_nb_vert_of_procs[ rank+1 ] += all_nb_elems[ i ];
}
- _nb_vert_of_procs[0] = 1; // base = 1
- for ( int i = 1; i < _nb_vert_of_procs.size(); ++i )
- _nb_vert_of_procs[ i ] += _nb_vert_of_procs[ i-1 ]; // to CSR format
- }
-// else
-// {
-// // to compute global ids of faces in joints
-// //_total_nb_faces = total_nb;
-// }
-
-// if ( !_rank) {
-// MEDMEM::STRING out("_nb_vert_of_procs :");
-// for ( int i = 0; i < _nb_vert_of_procs.size(); ++i )
-// out << _nb_vert_of_procs[i] << " ";
-// std::cout << out << std::endl;
-// }
-// if ( !_rank) {
-// MEDMEM::STRING out("elem_shift_by_domain :");
-// for ( int i = 0; i < elem_shift_by_domain.size(); ++i )
-// out << elem_shift_by_domain[i] << " ";
-// std::cout << out << std::endl;
-// }
-
+ _nb_vert_of_procs[0] = 0; // base = 0
+ for ( int i = 1; i < _nb_vert_of_procs.size(); ++i )
+ _nb_vert_of_procs[ i ] += _nb_vert_of_procs[ i-1 ]; // to CSR format
+
evaluateMemory();
-
+
return total_nb;
}
// 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();
// {
// vector<int> 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
const vector<int>& 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;
// 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();
// }
//================================================================================
}
return _max_memory - _init_memory;
}
+
+
+void ParaDomainSelector::sendMesh(const ParaMEDMEM::MEDCouplingUMesh& mesh, int target) const
+{
+
+
+ // First stage : sending sizes
+ // ------------------------------
+ vector<int> tinyInfoLocal;
+ vector<string> 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<int> tinyInfoDistant;
+ vector<string> 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<std::string> unusedTinyDistantSts;
+ mesh=dynamic_cast<ParaMEDMEM::MEDCouplingUMesh*> (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<ParaMEDMEM::MEDCouplingUMesh*> (distant_mesh_tmp);
+ //finish unserialization
+ mesh->unserialization(tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts);
+ std::cout<<"mesh size on recv"<<mesh->getNumberOfCells()<<std::endl;
+ //
+ if(v1Distant)
+ v1Distant->decrRef();
+ if(v2Distant)
+ v2Distant->decrRef();
+
+}
+void ParaDomainSelector::sendDoubleVec(const std::vector<double>& 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<double*>(&vec[0]), size,MPI_DOUBLE, target, 1212, MPI_COMM_WORLD);
+#endif
+}
+void ParaDomainSelector::recvDoubleVec(std::vector<double>& 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<int>& 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<int*>(&vec[0]), size,MPI_INT, target, 1212, MPI_COMM_WORLD);
+#endif
+}
+void ParaDomainSelector::recvIntVec(std::vector<int>& 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
+}
+
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())
// 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<double>& vec, int target) const;
+
+ void recvDoubleVec(std::vector<double>& vec, int source) const;
+
+ void sendIntVec(const std::vector<int>& vec, int target) const;
+ void recvIntVec(std::vector<int>& vec, int source) const;
+
private:
int _rank, _world_size; // my rank and nb of processors
MEDPARTITIONER_MESHCollectionMedAsciiDriver.hxx \
MEDPARTITIONER_ParallelTopology.hxx \
MEDPARTITIONER_FaceModel.hxx \
+MEDPARTITIONER_JointFinder.hxx \
MEDPARTITIONER_Graph.hxx\
MEDPARTITIONER_UserGraph.hxx\
MEDPARTITIONER_SequentialTopology.hxx \
MEDPARTITIONER_UserGraph.cxx\
MEDPARTITIONER_ParaDomainSelector.cxx \
MEDPARTITIONER_JointExchangeData.cxx \
+MEDPARTITIONER_JointFinder.cxx \
MEDPARTITIONER_SkyLineArray.cxx \
MEDPARTITIONER_ConnectZone.cxx
//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<std::string> unusedTinyDistantSts;
distant_mesh_tmp->resizeForUnserialization(tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts);
int nbLocalElems=0;
//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; i<triofield._nb_elems;i++)
{
for(int j=0;j<triofield._nodes_per_elem;j++)
}
delete[] conn;
- _mesh->setMeshDimension(triofield._mesh_dim);
_mesh->finishInsertingCells();
//field on the sending end
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) \
