--- /dev/null
+# Copyright (C) 2007-2017 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, or (at your option) any later version.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+# Author : Anthony Geay (EDF R&D)
+
+import re
+
+s1=1080
+s2=1104
+f=file("InterpKernelGaussCoords.cxx","r")
+lines=[elt[:-1] for elt in f.readlines()[s1:s2]]
+pat0=re.compile("void[\s]+GaussInfo\:\:([^\(]+)\([\s]*\)[\s]*$")
+pat1=re.compile("[\s]*\{[\s]*$")
+pat2=re.compile("[\s]+LOCAL_COORD_MACRO_BEGIN[\s]*\;[\s]*$")
+m0=pat0.match(lines[0])
+m1=pat1.match(lines[1])
+m2=pat2.match(lines[2])
+if (not m0) or (not m1) or (not m2):
+ raise Exception("Invalid first lines")
+offsetLines=3
+patEnd=re.compile("[\s]+LOCAL_COORD_MACRO_END[\s]*\;[\s]*$")
+mEnd=patEnd.match(lines[-1])
+if not mEnd:
+ raise Exception("Invalid end lines")
+#
+nbLines=len(lines)-4
+casePat=re.compile("[\s]+case[\s]+([\d]+)\:[\s]*$")
+entries=filter(lambda (i,x): casePat.match(x),enumerate(lines[offsetLines:-1]))
+#
+nbPts=len(entries)
+if nbLines%nbPts!=0:
+ raise Exception("Invalid lines nb !")
+dim=nbLines/nbPts-2
+if dim<1 or dim>3:
+ raise Exception("Ooops invalid dim !")
+entries=[(i,int(casePat.match(elt).group(1))) for i,elt in entries]
+assert(set([elt[1] for elt in entries])==set(range(nbPts)))
+#
+partEndEntries=re.compile("[\s]*break[\s]*\;[\s]*$")
+zePat=re.compile("[\s]+coords\[([\d]+)\][\s]*=[\s]*([\d]+[\.]?[\d]*)[\s]*\;[\s]*$")
+zeTab=(nbPts*dim)*[None]
+for lineId,ptId in entries:
+ endLine=lines[offsetLines+lineId+1+dim]
+ assert(partEndEntries.match(endLine))
+ for j in xrange(dim):
+ curLine=lines[offsetLines+lineId+1+j]
+ m=zePat.match(curLine)
+ assert(m)
+ assert(int(m.group(1))==j)
+ zeTab[ptId*dim+j]=m.group(2)
+ pass
+ pass
+assert(None not in zeTab)
+patInit="Init"
+assert(m0.group(1)[-len(patInit):]==patInit)
+varName="%s_REF"%((m0.group(1)[:-len(patInit)]).upper())
+print("const double %s[%d]={%s};"%(varName,len(zeTab),", ".join(zeTab)))
+for i in xrange(nbPts):
+ print(" case %d:"%(i))
+ for j in xrange(dim):
+ print(" coords[%d] = %s[%d];"%(j,varName,i*dim+j))
+ pass
+ print(" break;")
+
using namespace INTERP_KERNEL;
+const double GaussInfo::TETRA4A_REF[12]={0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0};
+
+const double GaussInfo::TETRA4B_REF[12]={0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0};
+
+const double GaussInfo::TETRA10A_REF[30]={0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.5, 0.5, 0.0, 0.5, 0.0, 0.5, 0.5, 0.0, 0.0};
+
+const double GaussInfo::TETRA10B_REF[30]={0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.5, 0.5, 0.5, 0.0, 0.5, 0.0, 0.0, 0.5, 0.0, 0.5};
+
+
//Define common part of the code in the MACRO
//---------------------------------------------------------------
#define LOCAL_COORD_MACRO_BEGIN \
return _my_nb_ref;
}
+GaussInfo GaussInfo::convertToLinear() const
+{
+ switch(_my_geometry)
+ {
+ case NORM_TETRA10:
+ {
+ std::vector<double> a(TETRA10A_REF,TETRA10A_REF+30),b(TETRA10B_REF,TETRA10B_REF+30);
+ if(IsSatisfy(a,_my_reference_coord))
+ {
+ std::vector<double> c(TETRA4A_REF,TETRA4A_REF+12);
+ return GaussInfo(NORM_TETRA4,_my_gauss_coord,getNbGauss(),c,4);
+ }
+ if(IsSatisfy(b,_my_reference_coord))
+ {
+ std::vector<double> c(TETRA4B_REF,TETRA4B_REF+12);
+ return GaussInfo(NORM_TETRA4,_my_gauss_coord,getNbGauss(),c,4);
+ }
+ throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for TETRA10 !");
+ }
+ default:
+ throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not implemented yet for other types than TETRA10 !");
+ }
+}
+
+
+bool GaussInfo::IsSatisfy(const std::vector<double>& ref1, const std::vector<double>& ref2)
+{
+ std::size_t sz(ref1.size());
+ if(sz!=ref2.size())
+ return false;
+ for(std::size_t i=0;i<sz;i++)
+ if(!IsEqual(ref1[i],ref2[i]))
+ return false;
+ return true;
+}
+
/*!
* Check coordinates
*/
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;
+ case 0:
+ coords[0] = TETRA4A_REF[0];
+ coords[1] = TETRA4A_REF[1];
+ coords[2] = TETRA4A_REF[2];
+ break;
+ case 1:
+ coords[0] = TETRA4A_REF[3];
+ coords[1] = TETRA4A_REF[4];
+ coords[2] = TETRA4A_REF[5];
+ break;
+ case 2:
+ coords[0] = TETRA4A_REF[6];
+ coords[1] = TETRA4A_REF[7];
+ coords[2] = TETRA4A_REF[8];
+ break;
+ case 3:
+ coords[0] = TETRA4A_REF[9];
+ coords[1] = TETRA4A_REF[10];
+ coords[2] = TETRA4A_REF[11];
+ 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;
}
/*!
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;
-
+ case 0:
+ coords[0] = TETRA4B_REF[0];
+ coords[1] = TETRA4B_REF[1];
+ coords[2] = TETRA4B_REF[2];
+ break;
+ case 1:
+ coords[0] = TETRA4B_REF[3];
+ coords[1] = TETRA4B_REF[4];
+ coords[2] = TETRA4B_REF[5];
+ break;
+ case 2:
+ coords[0] = TETRA4B_REF[6];
+ coords[1] = TETRA4B_REF[7];
+ coords[2] = TETRA4B_REF[8];
+ break;
+ case 3:
+ coords[0] = TETRA4B_REF[9];
+ coords[1] = TETRA4B_REF[10];
+ coords[2] = TETRA4B_REF[11];
+ 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;
}
/*!
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;
+ case 0:
+ coords[0] = TETRA10A_REF[0];
+ coords[1] = TETRA10A_REF[1];
+ coords[2] = TETRA10A_REF[2];
+ break;
+ case 1:
+ coords[0] = TETRA10A_REF[3];
+ coords[1] = TETRA10A_REF[4];
+ coords[2] = TETRA10A_REF[5];
+ break;
+ case 2:
+ coords[0] = TETRA10A_REF[6];
+ coords[1] = TETRA10A_REF[7];
+ coords[2] = TETRA10A_REF[8];
+ break;
+ case 3:
+ coords[0] = TETRA10A_REF[9];
+ coords[1] = TETRA10A_REF[10];
+ coords[2] = TETRA10A_REF[11];
+ break;
+ case 4:
+ coords[0] = TETRA10A_REF[12];
+ coords[1] = TETRA10A_REF[13];
+ coords[2] = TETRA10A_REF[14];
+ break;
+ case 5:
+ coords[0] = TETRA10A_REF[15];
+ coords[1] = TETRA10A_REF[16];
+ coords[2] = TETRA10A_REF[17];
+ break;
+ case 6:
+ coords[0] = TETRA10A_REF[18];
+ coords[1] = TETRA10A_REF[19];
+ coords[2] = TETRA10A_REF[20];
+ break;
+ case 7:
+ coords[0] = TETRA10A_REF[21];
+ coords[1] = TETRA10A_REF[22];
+ coords[2] = TETRA10A_REF[23];
+ break;
+ case 8:
+ coords[0] = TETRA10A_REF[24];
+ coords[1] = TETRA10A_REF[25];
+ coords[2] = TETRA10A_REF[26];
+ break;
+ case 9:
+ coords[0] = TETRA10A_REF[27];
+ coords[1] = TETRA10A_REF[28];
+ coords[2] = TETRA10A_REF[29];
+ 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;
+ 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;
}
/*!
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;
+ case 0:
+ coords[0] = TETRA10B_REF[0];
+ coords[1] = TETRA10B_REF[1];
+ coords[2] = TETRA10B_REF[2];
+ break;
+ case 1:
+ coords[0] = TETRA10B_REF[3];
+ coords[1] = TETRA10B_REF[4];
+ coords[2] = TETRA10B_REF[5];
+ break;
+ case 2:
+ coords[0] = TETRA10B_REF[6];
+ coords[1] = TETRA10B_REF[7];
+ coords[2] = TETRA10B_REF[8];
+ break;
+ case 3:
+ coords[0] = TETRA10B_REF[9];
+ coords[1] = TETRA10B_REF[10];
+ coords[2] = TETRA10B_REF[11];
+ break;
+ case 4:
+ coords[0] = TETRA10B_REF[12];
+ coords[1] = TETRA10B_REF[13];
+ coords[2] = TETRA10B_REF[14];
+ break;
+ case 5:
+ coords[0] = TETRA10B_REF[15];
+ coords[1] = TETRA10B_REF[16];
+ coords[2] = TETRA10B_REF[17];
+ break;
+ case 6:
+ coords[0] = TETRA10B_REF[18];
+ coords[1] = TETRA10B_REF[19];
+ coords[2] = TETRA10B_REF[20];
+ break;
+ case 7:
+ coords[0] = TETRA10B_REF[21];
+ coords[1] = TETRA10B_REF[22];
+ coords[2] = TETRA10B_REF[23];
+ break;
+ case 8:
+ coords[0] = TETRA10B_REF[24];
+ coords[1] = TETRA10B_REF[25];
+ coords[2] = TETRA10B_REF[26];
+ break;
+ case 9:
+ coords[0] = TETRA10B_REF[27];
+ coords[1] = TETRA10B_REF[28];
+ coords[2] = TETRA10B_REF[29];
+ 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;
}
/*!
INTERPKERNEL_EXPORT int getGaussCoordDim() const;
INTERPKERNEL_EXPORT int getReferenceCoordDim() const;
+ INTERPKERNEL_EXPORT DataVector getGaussCoords() const { return _my_gauss_coord; }
+ INTERPKERNEL_EXPORT DataVector getRefCoords() const { return _my_reference_coord; }
+ INTERPKERNEL_EXPORT NormalizedCellType getGeoType() const { return _my_geometry; }
INTERPKERNEL_EXPORT int getNbGauss() const;
INTERPKERNEL_EXPORT int getNbRef() const;
+ INTERPKERNEL_EXPORT GaussInfo convertToLinear() const;
+
INTERPKERNEL_EXPORT const double* getFunctionValues( const int theGaussId ) const;
INTERPKERNEL_EXPORT void initLocalInfo();
INTERPKERNEL_EXPORT static std::vector<double> NormalizeCoordinatesIfNecessary(NormalizedCellType ct, int inputDim, const std::vector<double>& inputArray);
- protected:
+ public:
+ static const double TETRA4A_REF[12];
+ static const double TETRA4B_REF[12];
+ static const double TETRA10A_REF[30];
+ static const double TETRA10B_REF[30];
+ protected:
+ static bool IsSatisfy(const std::vector<double>& ref1, const std::vector<double>& ref2);
bool isSatisfy();
void point1Init();
#include "MEDCouplingNatureOfField.hxx"
#include "InterpKernelAutoPtr.hxx"
+#include "InterpKernelGaussCoords.hxx"
#include <sstream>
#include <limits>
ret->copyAllTinyAttrFrom(this);
return ret;
}
+ case ON_GAUSS_PT:
+ {
+ const MEDCouplingMesh *mesh(getMesh());
+ if(!mesh)
+ throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::convertQuadraticCellsToLinear : null mesh !");
+ MCAuto<MEDCouplingUMesh> umesh(mesh->buildUnstructured());
+ std::set<INTERP_KERNEL::NormalizedCellType> gt(umesh->getAllGeoTypes());
+ MCAuto<MEDCouplingFieldDouble> ret(MEDCouplingFieldDouble::New(ON_GAUSS_PT));
+ //
+ const MEDCouplingFieldDiscretization *disc(getDiscretization());
+ const MEDCouplingFieldDiscretizationGauss *disc2(dynamic_cast<const MEDCouplingFieldDiscretizationGauss *>(disc));
+ if(!disc2)
+ throw INTERP_KERNEL::Exception("convertQuadraticCellsToLinear : Not a ON_GAUSS_PT field");
+ std::set<INTERP_KERNEL::NormalizedCellType> gt2(umesh->getAllGeoTypes());
+ const DataArrayDouble *arr(getArray());
+ std::vector< MCAuto<DataArrayInt> > cellIdsV;
+ std::vector< MCAuto<MEDCouplingUMesh> > meshesV;
+ std::vector< MEDCouplingGaussLocalization > glV;
+ bool isZipReq(false);
+ for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator it=gt.begin();it!=gt.end();it++)
+ {
+ const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel(*it));
+ MCAuto<DataArrayInt> cellIds(umesh->giveCellsWithType(*it));
+ cellIdsV.push_back(cellIds);
+ MCAuto<MEDCouplingUMesh> part(umesh->buildPartOfMySelf(cellIds->begin(),cellIds->end()));
+ int id(disc2->getGaussLocalizationIdOfOneType(*it));
+ const MEDCouplingGaussLocalization& gl(disc2->getGaussLocalization(id));
+ if(!cm.isQuadratic())
+ {
+ glV.push_back(gl);
+ }
+ else
+ {
+ isZipReq=true;
+ part->convertQuadraticCellsToLinear();
+ INTERP_KERNEL::GaussInfo gi(*it,gl.getGaussCoords(),gl.getNumberOfGaussPt(),gl.getRefCoords(),gl.getNumberOfPtsInRefCell());
+ INTERP_KERNEL::GaussInfo gi2(gi.convertToLinear());
+ MEDCouplingGaussLocalization gl2(gi2.getGeoType(),gi2.getRefCoords(),gi2.getGaussCoords(),gl.getWeights());
+ glV.push_back(gl2);
+ }
+ meshesV.push_back(part);
+ }
+ //
+ {
+ std::vector< const MEDCouplingUMesh * > meshesPtr(VecAutoToVecOfCstPt(meshesV));
+ umesh=MEDCouplingUMesh::MergeUMeshesOnSameCoords(meshesPtr);
+ std::vector< const DataArrayInt * > zeCellIds(VecAutoToVecOfCstPt(cellIdsV));
+ MCAuto<DataArrayInt> zeIds(DataArrayInt::Aggregate(zeCellIds));
+ umesh->renumberCells(zeIds->begin());
+ umesh->setName(mesh->getName());
+ }
+ //
+ if(isZipReq)
+ umesh->zipCoords();
+ ret->setArray(const_cast<DataArrayDouble *>(getArray()));
+ ret->setMesh(umesh);
+ for(std::vector< MEDCouplingGaussLocalization >::const_iterator it=glV.begin();it!=glV.end();it++)
+ ret->setGaussLocalizationOnType((*it).getType(),(*it).getRefCoords(),(*it).getGaussCoords(),(*it).getWeights());
+ ret->copyAllTinyAttrFrom(this);
+ ret->checkConsistencyLight();
+ return ret;
+ }
default:
- throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::convertQuadraticCellsToLinear : Only available for fields on nodes !");
+ throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::convertQuadraticCellsToLinear : Only available for fields on nodes and on cells !");
}
}
checkConsistencyLight();
if(!vor)
throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::voronoizeGen : null pointer !");
- const MEDCouplingMesh *inpMesh(getMesh());
+ MCAuto<MEDCouplingFieldDouble> fieldToWO;
+ const MEDCouplingMesh *inpMeshBase(getMesh());
+ MCAuto<MEDCouplingUMesh> inpMesh(inpMeshBase->buildUnstructured());
+ std::string meshName(inpMesh->getName());
+ if(!inpMesh->isPresenceOfQuadratic())
+ fieldToWO=clone(false);
+ else
+ {
+ fieldToWO=convertQuadraticCellsToLinear();
+ inpMeshBase=fieldToWO->getMesh();
+ inpMesh=inpMeshBase->buildUnstructured();
+ }
int nbCells(inpMesh->getNumberOfCells());
- const MEDCouplingFieldDiscretization *disc(getDiscretization());
+ const MEDCouplingFieldDiscretization *disc(fieldToWO->getDiscretization());
const MEDCouplingFieldDiscretizationGauss *disc2(dynamic_cast<const MEDCouplingFieldDiscretizationGauss *>(disc));
if(!disc2)
throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::voronoize2D : Not a ON_GAUSS_PT field");
MCAuto<DataArrayDouble> ptsInReal;
disc2->getCellIdsHavingGaussLocalization(i,ids);
{
- MCAuto<MEDCouplingUMesh> tmp4(inpMesh->buildUnstructured());
- MCAuto<MEDCouplingUMesh> subMesh(tmp4->buildPartOfMySelf(&ids[0],&ids[0]+ids.size()));
+ MCAuto<MEDCouplingUMesh> subMesh(inpMesh->buildPartOfMySelf(&ids[0],&ids[0]+ids.size()));
ptsInReal=gl.localizePtsInRefCooForEachCell(vorCellsForCurDisc->getCoords(),subMesh);
}
int nbPtsPerCell(vorCellsForCurDisc->getNumberOfNodes());
}
std::vector< const MEDCouplingUMesh * > cellsPtr(VecAutoToVecOfCstPt(cells));
MCAuto<MEDCouplingUMesh> outMesh(MEDCouplingUMesh::MergeUMeshes(cellsPtr));
+ outMesh->setName(meshName);
MCAuto<MEDCouplingFieldDouble> onCells(MEDCouplingFieldDouble::New(ON_CELLS));
onCells->setMesh(outMesh);
{
- MCAuto<DataArrayDouble> arr(getArray()->deepCopy());
+ MCAuto<DataArrayDouble> arr(fieldToWO->getArray()->deepCopy());
onCells->setArray(arr);
}
onCells->setTimeUnit(getTimeUnit());
std::vector<int> cut3DCurve(mDesc1->getNumberOfCells(),-2);
for(const int *it=cellIds1D->begin();it!=cellIds1D->end();it++)
cut3DCurve[*it]=-1;
- mDesc1->split3DCurveWithPlane(origin,vec,eps,cut3DCurve);
+ bool sameNbNodes;
+ {
+ int oldNbNodes(mDesc1->getNumberOfNodes());
+ mDesc1->split3DCurveWithPlane(origin,vec,eps,cut3DCurve);
+ sameNbNodes=(mDesc1->getNumberOfNodes()==oldNbNodes);
+ }
std::vector< std::pair<int,int> > cut3DSurf(mDesc2->getNumberOfCells());
AssemblyForSplitFrom3DCurve(cut3DCurve,nodes,mDesc2->getNodalConnectivity()->begin(),mDesc2->getNodalConnectivityIndex()->begin(),
mDesc1->getNodalConnectivity()->begin(),mDesc1->getNodalConnectivityIndex()->begin(),
std::vector<std::vector<int> > res;
buildSubCellsFromCut(cut3DSurf,desc2->begin(),descIndx2->begin(),mDesc1->getCoords()->begin(),eps,res);
std::size_t sz(res.size());
- if(res.size()==mDesc1->getNumberOfCells())
+ if(res.size()==mDesc1->getNumberOfCells() && sameNbNodes)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::clipSingle3DCellByPlane : cell is not clipped !");
for(std::size_t i=0;i<sz;i++)
{
conn2I->pushBackSilent(conn2->getNumberOfTuples());
ret2->setConnectivity(conn2,conn2I,true);
ret2->checkConsistencyLight();
- ret2->writeVTK("ret2.vtu");
ret2->orientCorrectlyPolyhedrons();
return ret2;
}
#
fieldOnCell=field.voronoize(1e-12) # hot point
fieldOnCell.checkConsistencyLight()
+ fieldOnCell.writeVTK("tt.vtu")
self.assertEqual(fieldOnCell.getMesh().getNumberOfCells(),7)
- self.assertEqual(fieldOnCell.getMesh().getNumberOfNodes(),32)
+ self.assertEqual(fieldOnCell.getMesh().getNumberOfNodes(),34)
self.assertTrue(fieldOnCell.getArray().isEqual(field.getArray(),1e-12))
- meaRef=DataArrayDouble([1.3,1.0,1.1975,1.36,1.4775,0.8,0.865])
+ meaRef=DataArrayDouble([1.3,0.9421428571428572,1.1975,1.36,1.4775,0.8,0.922857142857143])
mea=fieldOnCell.getMesh().getMeasureField(True).getArray()
self.assertTrue(mea.isEqual(meaRef,1e-12))# the first important test is here
self.assertEqual(field.getDiscretization().getNbOfGaussLocalization(),1)
self.assertTrue(a.isIota(7))# the second important test is here ! Check that Gauss points are inside the associated cell in fieldOnCell !
self.assertTrue(b.isIota(8))
#
- self.assertEqual(fieldOnCell.getMesh().buildDescendingConnectivity()[0].getNumberOfCells(),2*7+22)# last little test to reduce chance of errors. For humans there 21 but last tiny edge is split into 2 subedges due to alg
+ self.assertEqual(fieldOnCell.getMesh().buildDescendingConnectivity()[0].getNumberOfCells(),2*7+21)
pass
def testVoronoi3DSurf_1(self):
self.assertEqual(m2.getName(),"Mesh")
pass
+ def testVoronoi3D_3(self):
+ """Non regression test to check MEDCouplingUMesh::clipSingle3DCellByPlane"""
+ coo=DataArrayDouble([0.,1.,0.,0.,0.,0.,0.,0.,1.,1.,0.,0.],4,3)
+ m=MEDCouplingUMesh("mesh",3)
+ m.setCoords(coo) ; m.allocateCells()
+ m.insertNextCell(NORM_TETRA4,[0,2,3,1])
+ f=MEDCouplingFieldDouble(ON_GAUSS_PT)
+ f.setMesh(m) ; f.setName("field")
+ f.setGaussLocalizationOnType(NORM_TETRA4,[0.,1.,0.,0.,0.,0.,0.,0.,1.,1.,0.,0.],[0.1381966011250105, 0.1381966011250105, 0.1381966011250105, 0.1381966011250105, 0.1381966011250105, 0.5854101966249685, 0.1381966011250105, 0.5854101966249685, 0.1381966011250105, 0.5854101966249685, 0.1381966011250105, 0.1381966011250105], [0.041667,0.041667,0.041667,0.041667])
+ f.setArray(DataArrayDouble([0,1,2,3]))
+ f3=f.voronoize(1e-12)
+ ref=DataArrayDouble([0.047256836610416179,0.03980327668541684,0.039803276685416833,0.039803276685416833])
+ self.assertTrue(f3.getMesh().getMeasureField(False).getArray().isEqual(ref,1e-12))
+ self.assertTrue(f3.getArray().isEqual(DataArrayDouble([0,1,2,3]),1e-12))
+ pass
+
+ def testVoronoi3D_4(self):
+ """Idem testVoronoi3D_3 except that here quadratic cells are considered"""
+ coo=DataArrayDouble([0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,1.0,0.0,0.0,0.0,0.5,0.0,0.0,0.0,0.5,0.0,0.5,0.5,0.5,0.5,0.0,0.5,0.0,0.0,0.5,0.0,0.5],10,3)
+ m=MEDCouplingUMesh("mesh",3)
+ m.setCoords(coo) ; m.allocateCells()
+ m.insertNextCell(NORM_TETRA10,[0,1,2,3,4,5,6,7,8,9])
+ f=MEDCouplingFieldDouble(ON_GAUSS_PT)
+ f.setMesh(m) ; f.setName("field")
+ f.setGaussLocalizationOnType(NORM_TETRA10,[0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,1.0,0.0,0.0,0.0,0.5,0.0,0.0,0.0,0.5,0.0,0.5,0.5,0.5,0.5,0.0,0.5,0.0,0.0,0.5,0.0,0.5],[0.1381966011250105, 0.1381966011250105, 0.1381966011250105, 0.1381966011250105, 0.1381966011250105, 0.5854101966249685, 0.1381966011250105, 0.5854101966249685, 0.1381966011250105, 0.5854101966249685, 0.1381966011250105, 0.1381966011250105], [0.041667,0.041667,0.041667,0.041667])
+ f.setArray(DataArrayDouble([0,1,2,3]))
+ f3=f.voronoize(1e-12)
+ ref=DataArrayDouble([0.047256836610416179,0.03980327668541684,0.039803276685416833,0.039803276685416833])
+ self.assertTrue(f3.getMesh().getMeasureField(False).getArray().isEqual(ref,1e-12))
+ self.assertTrue(f3.getArray().isEqual(DataArrayDouble([0,1,2,3]),1e-12))
+ pass
+
+ def testConvertQuadToLin4Gauss_1(self):
+ coo=DataArrayDouble([0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,1.0,0.0,0.0,0.0,0.5,0.0,0.0,0.0,0.5,0.0,0.5,0.5,0.5,0.5,0.0,0.5,0.0,0.0,0.5,0.0,0.5],10,3)
+ m=MEDCouplingUMesh("mesh",3)
+ m.setCoords(coo) ; m.allocateCells()
+ m.insertNextCell(NORM_TETRA10,[0,1,2,3,4,5,6,7,8,9])
+ f=MEDCouplingFieldDouble(ON_GAUSS_PT)
+ f.setMesh(m) ; f.setName("field")
+ aaaa=[0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,1.0,0.0,0.0,0.0,0.5,0.0,0.0,0.0,0.5,0.0,0.5,0.5,0.5,0.5,0.0,0.5,0.0,0.0,0.5,0.0,0.5]
+ bbbb=[0.1381966011250105,0.1381966011250105,0.1381966011250105,0.1381966011250105,0.1381966011250105,0.5854101966249685,0.1381966011250105,0.5854101966249685,0.1381966011250105,0.5854101966249685,0.1381966011250105,0.1381966011250105]
+ cccc=[0.041667,0.041667,0.041667,0.041667]
+ f.setGaussLocalizationOnType(NORM_TETRA10,aaaa,bbbb,cccc)
+ f.setArray(DataArrayDouble([0,1,2,3]))
+ f.setTime(1.,2,3)
+ #
+ mcpy=m.deepCopy() ; mcpy.convertQuadraticCellsToLinear() ; mcpy.zipCoords()
+ #
+ f2=f.convertQuadraticCellsToLinear()
+ f2.checkConsistencyLight()
+ self.assertTrue(f2.getMesh().isEqual(mcpy,1e-12))
+ self.assertTrue(f2.getArray().isEqual(DataArrayDouble([0,1,2,3]),1e-12))
+ self.assertEqual(f2.getNbOfGaussLocalization(),1)
+ gl=f2.getGaussLocalization(0)
+ self.assertEqual(gl.getType(),NORM_TETRA4)
+ self.assertTrue(DataArrayDouble(gl.getRefCoords()).isEqual(DataArrayDouble([0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,1.0,0.0,0.0]),1e-12))
+ self.assertTrue(DataArrayDouble(gl.getGaussCoords()).isEqual(DataArrayDouble(bbbb),1e-12))
+ self.assertTrue(DataArrayDouble(gl.getWeights()).isEqual(DataArrayDouble(cccc),1e-12))
+ self.assertEqual(f2.getName(),"field")
+ self.assertEqual(f2.getTime(),[1.,2,3])
+ pass
+
pass
if __name__ == '__main__':
