1 // Copyright (C) 2007-2012 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (CEA/DEN)
21 #include "MEDCouplingFieldDiscretization.hxx"
22 #include "MEDCouplingCMesh.hxx"
23 #include "MEDCouplingUMesh.hxx"
24 #include "MEDCouplingFieldDouble.hxx"
25 #include "MEDCouplingAutoRefCountObjectPtr.hxx"
27 #include "CellModel.hxx"
28 #include "InterpolationUtils.hxx"
29 #include "InterpKernelAutoPtr.hxx"
30 #include "InterpKernelGaussCoords.hxx"
31 #include "InterpKernelMatrixTools.hxx"
40 using namespace ParaMEDMEM;
42 const double MEDCouplingFieldDiscretization::DFLT_PRECISION=1.e-12;
44 const char MEDCouplingFieldDiscretizationP0::REPR[]="P0";
46 const TypeOfField MEDCouplingFieldDiscretizationP0::TYPE=ON_CELLS;
48 const char MEDCouplingFieldDiscretizationP1::REPR[]="P1";
50 const TypeOfField MEDCouplingFieldDiscretizationP1::TYPE=ON_NODES;
52 const int MEDCouplingFieldDiscretizationPerCell::DFT_INVALID_LOCID_VALUE=-1;
54 const char MEDCouplingFieldDiscretizationGauss::REPR[]="GAUSS";
56 const TypeOfField MEDCouplingFieldDiscretizationGauss::TYPE=ON_GAUSS_PT;
58 const char MEDCouplingFieldDiscretizationGaussNE::REPR[]="GSSNE";
60 const TypeOfField MEDCouplingFieldDiscretizationGaussNE::TYPE=ON_GAUSS_NE;
62 const char MEDCouplingFieldDiscretizationKriging::REPR[]="KRIGING";
64 const TypeOfField MEDCouplingFieldDiscretizationKriging::TYPE=ON_NODES_KR;
66 MEDCouplingFieldDiscretization::MEDCouplingFieldDiscretization():_precision(DFLT_PRECISION)
70 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretization::New(TypeOfField type)
74 case MEDCouplingFieldDiscretizationP0::TYPE:
75 return new MEDCouplingFieldDiscretizationP0;
76 case MEDCouplingFieldDiscretizationP1::TYPE:
77 return new MEDCouplingFieldDiscretizationP1;
78 case MEDCouplingFieldDiscretizationGauss::TYPE:
79 return new MEDCouplingFieldDiscretizationGauss;
80 case MEDCouplingFieldDiscretizationGaussNE::TYPE:
81 return new MEDCouplingFieldDiscretizationGaussNE;
82 case MEDCouplingFieldDiscretizationKriging::TYPE:
83 return new MEDCouplingFieldDiscretizationKriging;
85 throw INTERP_KERNEL::Exception("Choosen discretization is not implemented yet.");
89 TypeOfField MEDCouplingFieldDiscretization::getTypeOfFieldFromStringRepr(const char *repr) throw(INTERP_KERNEL::Exception)
91 std::string reprCpp(repr);
92 if(reprCpp==MEDCouplingFieldDiscretizationP0::REPR)
93 return MEDCouplingFieldDiscretizationP0::TYPE;
94 if(reprCpp==MEDCouplingFieldDiscretizationP1::REPR)
95 return MEDCouplingFieldDiscretizationP1::TYPE;
96 if(reprCpp==MEDCouplingFieldDiscretizationGauss::REPR)
97 return MEDCouplingFieldDiscretizationGauss::TYPE;
98 if(reprCpp==MEDCouplingFieldDiscretizationGaussNE::REPR)
99 return MEDCouplingFieldDiscretizationGaussNE::TYPE;
100 if(reprCpp==MEDCouplingFieldDiscretizationKriging::REPR)
101 return MEDCouplingFieldDiscretizationKriging::TYPE;
102 throw INTERP_KERNEL::Exception("Representation does not match with any field discretization !");
105 bool MEDCouplingFieldDiscretization::isEqual(const MEDCouplingFieldDiscretization *other, double eps) const
108 return isEqualIfNotWhy(other,eps,reason);
111 bool MEDCouplingFieldDiscretization::isEqualWithoutConsideringStr(const MEDCouplingFieldDiscretization *other, double eps) const
113 return isEqual(other,eps);
117 * For all field discretization excepted GaussPts the [ \a startCellIds, \a endCellIds ) has no impact on the cloned instance.
119 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretization::clonePart(const int *startCellIds, const int *endCellIds) const
125 * Excepted for MEDCouplingFieldDiscretizationPerCell no underlying TimeLabel object : nothing to do in generally.
127 void MEDCouplingFieldDiscretization::updateTime() const
131 std::size_t MEDCouplingFieldDiscretization::getHeapMemorySize() const
137 * Computes normL1 of DataArrayDouble instance arr.
138 * @param res output parameter expected to be of size arr->getNumberOfComponents();
139 * @throw when the field discretization fails on getMeasure fields (gauss points for example)
141 void MEDCouplingFieldDiscretization::normL1(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, double *res) const throw(INTERP_KERNEL::Exception)
143 MEDCouplingFieldDouble *vol=getMeasureField(mesh,true);
144 int nbOfCompo=arr->getNumberOfComponents();
145 int nbOfElems=getNumberOfTuples(mesh);
146 std::fill(res,res+nbOfCompo,0.);
147 const double *arrPtr=arr->getConstPointer();
148 const double *volPtr=vol->getArray()->getConstPointer();
150 for(int i=0;i<nbOfElems;i++)
152 double v=fabs(volPtr[i]);
153 for(int j=0;j<nbOfCompo;j++)
154 res[j]+=fabs(arrPtr[i*nbOfCompo+j])*v;
157 std::transform(res,res+nbOfCompo,res,std::bind2nd(std::multiplies<double>(),1./deno));
162 * Computes normL2 of DataArrayDouble instance arr.
163 * @param res output parameter expected to be of size arr->getNumberOfComponents();
164 * @throw when the field discretization fails on getMeasure fields (gauss points for example)
166 void MEDCouplingFieldDiscretization::normL2(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, double *res) const throw(INTERP_KERNEL::Exception)
168 MEDCouplingFieldDouble *vol=getMeasureField(mesh,true);
169 int nbOfCompo=arr->getNumberOfComponents();
170 int nbOfElems=getNumberOfTuples(mesh);
171 std::fill(res,res+nbOfCompo,0.);
172 const double *arrPtr=arr->getConstPointer();
173 const double *volPtr=vol->getArray()->getConstPointer();
175 for(int i=0;i<nbOfElems;i++)
177 double v=fabs(volPtr[i]);
178 for(int j=0;j<nbOfCompo;j++)
179 res[j]+=arrPtr[i*nbOfCompo+j]*arrPtr[i*nbOfCompo+j]*v;
182 std::transform(res,res+nbOfCompo,res,std::bind2nd(std::multiplies<double>(),1./deno));
183 std::transform(res,res+nbOfCompo,res,std::ptr_fun<double,double>(std::sqrt));
188 * Computes integral of DataArrayDouble instance arr.
189 * @param res output parameter expected to be of size arr->getNumberOfComponents();
190 * @throw when the field discretization fails on getMeasure fields (gauss points for example)
192 void MEDCouplingFieldDiscretization::integral(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception)
194 MEDCouplingFieldDouble *vol=getMeasureField(mesh,isWAbs);
195 int nbOfCompo=arr->getNumberOfComponents();
196 int nbOfElems=getNumberOfTuples(mesh);
197 std::fill(res,res+nbOfCompo,0.);
198 const double *arrPtr=arr->getConstPointer();
199 const double *volPtr=vol->getArray()->getConstPointer();
200 double *tmp=new double[nbOfCompo];
201 for (int i=0;i<nbOfElems;i++)
203 std::transform(arrPtr+i*nbOfCompo,arrPtr+(i+1)*nbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),volPtr[i]));
204 std::transform(tmp,tmp+nbOfCompo,res,res,std::plus<double>());
210 void MEDCouplingFieldDiscretization::getSerializationIntArray(DataArrayInt *& arr) const
218 void MEDCouplingFieldDiscretization::getTinySerializationIntInformation(std::vector<int>& tinyInfo) const
225 void MEDCouplingFieldDiscretization::getTinySerializationDbleInformation(std::vector<double>& tinyInfo) const
229 void MEDCouplingFieldDiscretization::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *& arr)
237 void MEDCouplingFieldDiscretization::finishUnserialization(const std::vector<double>& tinyInfo)
242 * This method is typically the first step of renumbering. The implementation is empty it is not a bug only gauss is impacted
243 * virtualy by this method.
245 void MEDCouplingFieldDiscretization::renumberCells(const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
249 double MEDCouplingFieldDiscretization::getIJK(const MEDCouplingMesh *mesh, const DataArrayDouble *da,
250 int cellId, int nodeIdInCell, int compoId) const throw(INTERP_KERNEL::Exception)
252 throw INTERP_KERNEL::Exception("getIJK Invalid ! only for GaussPoint and GaussNE discretizations !");
255 void MEDCouplingFieldDiscretization::setGaussLocalizationOnType(const MEDCouplingMesh *m, INTERP_KERNEL::NormalizedCellType type, const std::vector<double>& refCoo,
256 const std::vector<double>& gsCoo, const std::vector<double>& wg) throw(INTERP_KERNEL::Exception)
258 throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
261 void MEDCouplingFieldDiscretization::setGaussLocalizationOnCells(const MEDCouplingMesh *m, const int *begin, const int *end, const std::vector<double>& refCoo,
262 const std::vector<double>& gsCoo, const std::vector<double>& wg) throw(INTERP_KERNEL::Exception)
264 throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
267 void MEDCouplingFieldDiscretization::clearGaussLocalizations() throw(INTERP_KERNEL::Exception)
269 throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
272 MEDCouplingGaussLocalization& MEDCouplingFieldDiscretization::getGaussLocalization(int locId) throw(INTERP_KERNEL::Exception)
274 throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
277 const MEDCouplingGaussLocalization& MEDCouplingFieldDiscretization::getGaussLocalization(int locId) const throw(INTERP_KERNEL::Exception)
279 throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
282 int MEDCouplingFieldDiscretization::getNbOfGaussLocalization() const throw(INTERP_KERNEL::Exception)
284 throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
287 int MEDCouplingFieldDiscretization::getGaussLocalizationIdOfOneCell(int cellId) const throw(INTERP_KERNEL::Exception)
289 throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
292 int MEDCouplingFieldDiscretization::getGaussLocalizationIdOfOneType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception)
294 throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
297 std::set<int> MEDCouplingFieldDiscretization::getGaussLocalizationIdsOfOneType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception)
299 throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
302 void MEDCouplingFieldDiscretization::getCellIdsHavingGaussLocalization(int locId, std::vector<int>& cellIds) const throw(INTERP_KERNEL::Exception)
304 throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
307 void MEDCouplingFieldDiscretization::RenumberEntitiesFromO2NArr(double eps, const int *old2NewPtr, int newNbOfEntity, DataArrayDouble *arr, const char *msg)
309 int oldNbOfElems=arr->getNumberOfTuples();
310 int nbOfComp=arr->getNumberOfComponents();
311 int newNbOfTuples=newNbOfEntity;
312 DataArrayDouble *arrCpy=arr->deepCpy();
313 const double *ptSrc=arrCpy->getConstPointer();
314 arr->reAlloc(newNbOfTuples);
315 double *ptToFill=arr->getPointer();
316 std::fill(ptToFill,ptToFill+nbOfComp*newNbOfTuples,std::numeric_limits<double>::max());
317 INTERP_KERNEL::AutoPtr<double> tmp=new double[nbOfComp];
318 for(int i=0;i<oldNbOfElems;i++)
320 int newNb=old2NewPtr[i];
321 if(newNb>=0)//if newNb<0 the node is considered as out.
323 if(std::find_if(ptToFill+newNb*nbOfComp,ptToFill+(newNb+1)*nbOfComp,std::bind2nd(std::not_equal_to<double>(),std::numeric_limits<double>::max()))
324 ==ptToFill+(newNb+1)*nbOfComp)
325 std::copy(ptSrc+i*nbOfComp,ptSrc+(i+1)*nbOfComp,ptToFill+newNb*nbOfComp);
328 std::transform(ptSrc+i*nbOfComp,ptSrc+(i+1)*nbOfComp,ptToFill+newNb*nbOfComp,(double *)tmp,std::minus<double>());
329 std::transform((double *)tmp,((double *)tmp)+nbOfComp,(double *)tmp,std::ptr_fun<double,double>(fabs));
330 //if(!std::equal(ptSrc+i*nbOfComp,ptSrc+(i+1)*nbOfComp,ptToFill+newNb*nbOfComp))
331 if(*std::max_element((double *)tmp,((double *)tmp)+nbOfComp)>eps)
334 std::ostringstream oss;
335 oss << msg << " " << i << " and " << std::find(old2NewPtr,old2NewPtr+i,newNb)-old2NewPtr
336 << " have been merged and " << msg << " field on them are different !";
337 throw INTERP_KERNEL::Exception(oss.str().c_str());
345 void MEDCouplingFieldDiscretization::RenumberEntitiesFromN2OArr(const int *new2OldPtr, int new2OldSz, DataArrayDouble *arr, const char *msg)
347 int nbOfComp=arr->getNumberOfComponents();
348 DataArrayDouble *arrCpy=arr->deepCpy();
349 const double *ptSrc=arrCpy->getConstPointer();
350 arr->reAlloc(new2OldSz);
351 double *ptToFill=arr->getPointer();
352 for(int i=0;i<new2OldSz;i++)
354 int oldNb=new2OldPtr[i];
355 std::copy(ptSrc+oldNb*nbOfComp,ptSrc+(oldNb+1)*nbOfComp,ptToFill+i*nbOfComp);
360 MEDCouplingFieldDiscretization::~MEDCouplingFieldDiscretization()
364 TypeOfField MEDCouplingFieldDiscretizationP0::getEnum() const
369 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationP0::clone() const
371 return new MEDCouplingFieldDiscretizationP0;
374 std::string MEDCouplingFieldDiscretizationP0::getStringRepr() const
376 return std::string(REPR);
379 const char *MEDCouplingFieldDiscretizationP0::getRepr() const
384 bool MEDCouplingFieldDiscretizationP0::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
386 const MEDCouplingFieldDiscretizationP0 *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationP0 *>(other);
389 reason="Spatial discrtization of this is ON_CELLS, which is not the case of other.";
393 int MEDCouplingFieldDiscretizationP0::getNumberOfTuples(const MEDCouplingMesh *mesh) const
395 return mesh->getNumberOfCells();
398 int MEDCouplingFieldDiscretizationP0::getNumberOfMeshPlaces(const MEDCouplingMesh *mesh) const
400 return mesh->getNumberOfCells();
403 DataArrayInt *MEDCouplingFieldDiscretizationP0::getOffsetArr(const MEDCouplingMesh *mesh) const
405 int nbOfTuples=mesh->getNumberOfCells();
406 DataArrayInt *ret=DataArrayInt::New();
407 ret->alloc(nbOfTuples+1,1);
412 void MEDCouplingFieldDiscretizationP0::renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArrayDouble *>& arrays,
413 const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
415 const int *array=old2NewBg;
417 array=DataArrayInt::CheckAndPreparePermutation(old2NewBg,old2NewBg+mesh->getNumberOfCells());
418 for(std::vector<DataArrayDouble *>::const_iterator it=arrays.begin();it!=arrays.end();it++)
421 (*it)->renumberInPlace(array);
427 DataArrayDouble *MEDCouplingFieldDiscretizationP0::getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const
429 return mesh->getBarycenterAndOwner();
432 void MEDCouplingFieldDiscretizationP0::computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const int *partBg, const int *partEnd,
433 DataArrayInt *&cellRest)
435 cellRest=DataArrayInt::New();
436 cellRest->alloc((int)std::distance(partBg,partEnd),1);
437 std::copy(partBg,partEnd,cellRest->getPointer());
440 void MEDCouplingFieldDiscretizationP0::checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception)
444 void MEDCouplingFieldDiscretizationP0::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception)
446 if(mesh->getNumberOfCells()!=da->getNumberOfTuples())
448 std::ostringstream message;
449 message << "Field on cells invalid because there are " << mesh->getNumberOfCells();
450 message << " cells in mesh and " << da->getNumberOfTuples() << " tuples in field !";
451 throw INTERP_KERNEL::Exception(message.str().c_str());
455 MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP0::getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const
458 throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP0::getMeasureField : mesh instance specified is NULL !");
459 return mesh->getMeasureField(isAbs);
462 void MEDCouplingFieldDiscretizationP0::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
464 int id=mesh->getCellContainingPoint(loc,_precision);
466 throw INTERP_KERNEL::Exception("Specified point is detected outside of mesh : unable to apply P0::getValueOn !");
467 arr->getTuple(id,res);
470 void MEDCouplingFieldDiscretizationP0::getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const
472 const MEDCouplingCMesh *meshC=dynamic_cast<const MEDCouplingCMesh *>(mesh);
474 throw INTERP_KERNEL::Exception("P0::getValueOnPos is only accessible for structured meshes !");
475 int id=meshC->getCellIdFromPos(i,j,k);
476 arr->getTuple(id,res);
479 DataArrayDouble *MEDCouplingFieldDiscretizationP0::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfPoints) const
481 std::vector<int> elts,eltsIndex;
482 mesh->getCellsContainingPoints(loc,nbOfPoints,_precision,elts,eltsIndex);
483 int spaceDim=mesh->getSpaceDimension();
484 int nbOfComponents=arr->getNumberOfComponents();
485 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=DataArrayDouble::New();
486 ret->alloc(nbOfPoints,nbOfComponents);
487 double *ptToFill=ret->getPointer();
488 for(int i=0;i<nbOfPoints;i++,ptToFill+=nbOfComponents)
489 if(eltsIndex[i+1]-eltsIndex[i]>=1)
490 arr->getTuple(elts[eltsIndex[i]],ptToFill);
493 std::ostringstream oss; oss << "Point #" << i << " with coordinates : (";
494 std::copy(loc+i*spaceDim,loc+(i+1)*spaceDim,std::ostream_iterator<double>(oss,", "));
495 oss << ") detected outside mesh : unable to apply P0::getValueOnMulti ! ";
496 throw INTERP_KERNEL::Exception(oss.str().c_str());
502 * Nothing to do. It's not a bug.
504 void MEDCouplingFieldDiscretizationP0::renumberValuesOnNodes(double , const int *, int newNbOfNodes, DataArrayDouble *) const
508 void MEDCouplingFieldDiscretizationP0::renumberValuesOnCells(double epsOnVals, const MEDCouplingMesh *mesh, const int *old2New, int newSz, DataArrayDouble *arr) const
510 RenumberEntitiesFromO2NArr(epsOnVals,old2New,newSz,arr,"Cell");
513 void MEDCouplingFieldDiscretizationP0::renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const
515 RenumberEntitiesFromN2OArr(new2old,newSz,arr,"Cell");
519 * This method returns a tuple ids selection from cell ids selection [start;end).
520 * This method is called by MEDCouplingFieldDiscretizationP0::buildSubMeshData to return parameter \b di.
521 * Here for P0 it's very simple !
523 * \return a newly allocated array containing ids to select into the DataArrayDouble of the field.
526 DataArrayInt *MEDCouplingFieldDiscretizationP0::computeTupleIdsToSelectFromCellIds(const MEDCouplingMesh *mesh, const int *startCellIds, const int *endCellIds) const
528 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
529 ret->alloc((int)std::distance(startCellIds,endCellIds),1);
530 std::copy(startCellIds,endCellIds,ret->getPointer());
535 * This method returns a submesh of 'mesh' instance constituting cell ids contained in array defined as an interval [start;end).
536 * @param di is an array returned that specifies entity ids (here cells ids) in mesh 'mesh' of entity in returned submesh.
537 * Example : The first cell id of returned mesh has the (*di)[0] id in 'mesh'
539 MEDCouplingMesh *MEDCouplingFieldDiscretizationP0::buildSubMeshData(const MEDCouplingMesh *mesh, const int *start, const int *end, DataArrayInt *&di) const
541 MEDCouplingMesh *ret=mesh->buildPart(start,end);
542 di=DataArrayInt::New();
543 di->alloc((int)std::distance(start,end),1);
544 int *pt=di->getPointer();
545 std::copy(start,end,pt);
549 int MEDCouplingFieldDiscretizationOnNodes::getNumberOfTuples(const MEDCouplingMesh *mesh) const
551 return mesh->getNumberOfNodes();
554 int MEDCouplingFieldDiscretizationOnNodes::getNumberOfMeshPlaces(const MEDCouplingMesh *mesh) const
556 return mesh->getNumberOfNodes();
560 * Nothing to do here.
562 void MEDCouplingFieldDiscretizationOnNodes::renumberArraysForCell(const MEDCouplingMesh *, const std::vector<DataArrayDouble *>& arrays,
563 const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
567 DataArrayInt *MEDCouplingFieldDiscretizationOnNodes::getOffsetArr(const MEDCouplingMesh *mesh) const
569 int nbOfTuples=mesh->getNumberOfNodes();
570 DataArrayInt *ret=DataArrayInt::New();
571 ret->alloc(nbOfTuples+1,1);
576 DataArrayDouble *MEDCouplingFieldDiscretizationOnNodes::getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const
578 return mesh->getCoordinatesAndOwner();
581 void MEDCouplingFieldDiscretizationOnNodes::computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const int *partBg, const int *partEnd,
582 DataArrayInt *&cellRest)
584 cellRest=mesh->getCellIdsFullyIncludedInNodeIds(partBg,partEnd);
587 void MEDCouplingFieldDiscretizationOnNodes::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception)
589 if(mesh->getNumberOfNodes()!=da->getNumberOfTuples())
591 std::ostringstream message;
592 message << "Field on nodes invalid because there are " << mesh->getNumberOfNodes();
593 message << " nodes in mesh and " << da->getNumberOfTuples() << " tuples in field !";
594 throw INTERP_KERNEL::Exception(message.str().c_str());
599 * This method returns a submesh of 'mesh' instance constituting cell ids contained in array defined as an interval [start;end).
600 * @param di is an array returned that specifies entity ids (here nodes ids) in mesh 'mesh' of entity in returned submesh.
601 * Example : The first node id of returned mesh has the (*di)[0] id in 'mesh'
603 MEDCouplingMesh *MEDCouplingFieldDiscretizationOnNodes::buildSubMeshData(const MEDCouplingMesh *mesh, const int *start, const int *end, DataArrayInt *&di) const
605 MEDCouplingMesh *ret=mesh->buildPartAndReduceNodes(start,end,di);
606 DataArrayInt *di2=di->invertArrayO2N2N2O(ret->getNumberOfNodes());
613 * This method returns a tuple ids selection from cell ids selection [start;end).
614 * This method is called by MEDCouplingFieldDiscretizationP0::buildSubMeshData to return parameter \b di.
615 * Here for P1 only nodes fetched by submesh of mesh[startCellIds:endCellIds) is returned !
617 * \return a newly allocated array containing ids to select into the DataArrayDouble of the field.
620 DataArrayInt *MEDCouplingFieldDiscretizationOnNodes::computeTupleIdsToSelectFromCellIds(const MEDCouplingMesh *mesh, const int *startCellIds, const int *endCellIds) const
623 throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP1::computeTupleIdsToSelectFromCellIds : null mesh !");
624 const MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> umesh=mesh->buildUnstructured();
625 MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> umesh2=static_cast<MEDCouplingUMesh *>(umesh->buildPartOfMySelf(startCellIds,endCellIds,true));
626 return umesh2->computeFetchedNodeIds();
629 void MEDCouplingFieldDiscretizationOnNodes::renumberValuesOnNodes(double epsOnVals, const int *old2NewPtr, int newNbOfNodes, DataArrayDouble *arr) const
631 RenumberEntitiesFromO2NArr(epsOnVals,old2NewPtr,newNbOfNodes,arr,"Node");
635 * Nothing to do it's not a bug.
637 void MEDCouplingFieldDiscretizationOnNodes::renumberValuesOnCells(double epsOnVals, const MEDCouplingMesh *mesh, const int *old2New, int newSz, DataArrayDouble *arr) const
642 * Nothing to do it's not a bug.
644 void MEDCouplingFieldDiscretizationOnNodes::renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const
648 void MEDCouplingFieldDiscretizationOnNodes::getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const
650 const MEDCouplingCMesh *meshC=dynamic_cast<const MEDCouplingCMesh *>(mesh);
652 throw INTERP_KERNEL::Exception("OnNodes::getValueOnPos(i,j,k) is only accessible for structured meshes !");
653 int id=meshC->getNodeIdFromPos(i,j,k);
654 arr->getTuple(id,res);
657 TypeOfField MEDCouplingFieldDiscretizationP1::getEnum() const
662 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationP1::clone() const
664 return new MEDCouplingFieldDiscretizationP1;
667 std::string MEDCouplingFieldDiscretizationP1::getStringRepr() const
669 return std::string(REPR);
672 const char *MEDCouplingFieldDiscretizationP1::getRepr() const
677 bool MEDCouplingFieldDiscretizationP1::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
679 const MEDCouplingFieldDiscretizationP1 *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationP1 *>(other);
682 reason="Spatial discrtization of this is ON_NODES, which is not the case of other.";
686 void MEDCouplingFieldDiscretizationP1::checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception)
688 if(nat!=ConservativeVolumic)
689 throw INTERP_KERNEL::Exception("Invalid nature for P1 field : expected ConservativeVolumic !");
692 MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP1::getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const
695 throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP1::getMeasureField : mesh instance specified is NULL !");
696 return mesh->getMeasureFieldOnNode(isAbs);
699 void MEDCouplingFieldDiscretizationP1::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
701 int id=mesh->getCellContainingPoint(loc,_precision);
703 throw INTERP_KERNEL::Exception("Specified point is detected outside of mesh : unable to apply P1::getValueOn !");
704 INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell(id);
705 if(type!=INTERP_KERNEL::NORM_SEG2 && type!=INTERP_KERNEL::NORM_TRI3 && type!=INTERP_KERNEL::NORM_TETRA4)
706 throw INTERP_KERNEL::Exception("P1 getValueOn is not specified for not simplex cells !");
707 getValueInCell(mesh,id,arr,loc,res);
711 * This method localizes a point defined by 'loc' in a cell with id 'cellId' into mesh 'mesh'.
712 * The result is put into res expected to be of size at least arr->getNumberOfComponents()
714 void MEDCouplingFieldDiscretizationP1::getValueInCell(const MEDCouplingMesh *mesh, int cellId, const DataArrayDouble *arr, const double *loc, double *res) const
716 std::vector<int> conn;
717 std::vector<double> coo;
718 mesh->getNodeIdsOfCell(cellId,conn);
719 for(std::vector<int>::const_iterator iter=conn.begin();iter!=conn.end();iter++)
720 mesh->getCoordinatesOfNode(*iter,coo);
721 int spaceDim=mesh->getSpaceDimension();
722 std::size_t nbOfNodes=conn.size();
723 std::vector<const double *> vec(nbOfNodes);
724 for(std::size_t i=0;i<nbOfNodes;i++)
725 vec[i]=&coo[i*spaceDim];
726 INTERP_KERNEL::AutoPtr<double> tmp=new double[nbOfNodes];
727 INTERP_KERNEL::barycentric_coords(vec,loc,tmp);
728 int sz=arr->getNumberOfComponents();
729 INTERP_KERNEL::AutoPtr<double> tmp2=new double[sz];
730 std::fill(res,res+sz,0.);
731 for(std::size_t i=0;i<nbOfNodes;i++)
733 arr->getTuple(conn[i],(double *)tmp2);
734 std::transform((double *)tmp2,((double *)tmp2)+sz,(double *)tmp2,std::bind2nd(std::multiplies<double>(),tmp[i]));
735 std::transform(res,res+sz,(double *)tmp2,res,std::plus<double>());
739 DataArrayDouble *MEDCouplingFieldDiscretizationP1::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfPoints) const
741 std::vector<int> elts,eltsIndex;
742 mesh->getCellsContainingPoints(loc,nbOfPoints,_precision,elts,eltsIndex);
743 int spaceDim=mesh->getSpaceDimension();
744 int nbOfComponents=arr->getNumberOfComponents();
745 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=DataArrayDouble::New();
746 ret->alloc(nbOfPoints,nbOfComponents);
747 double *ptToFill=ret->getPointer();
748 for(int i=0;i<nbOfPoints;i++)
749 if(eltsIndex[i+1]-eltsIndex[i]>=1)
750 getValueInCell(mesh,elts[eltsIndex[i]],arr,loc+i*spaceDim,ptToFill+i*nbOfComponents);
753 std::ostringstream oss; oss << "Point #" << i << " with coordinates : (";
754 std::copy(loc+i*spaceDim,loc+(i+1)*spaceDim,std::ostream_iterator<double>(oss,", "));
755 oss << ") detected outside mesh : unable to apply P1::getValueOnMulti ! ";
756 throw INTERP_KERNEL::Exception(oss.str().c_str());
761 MEDCouplingFieldDiscretizationPerCell::MEDCouplingFieldDiscretizationPerCell():_discr_per_cell(0)
765 MEDCouplingFieldDiscretizationPerCell::~MEDCouplingFieldDiscretizationPerCell()
768 _discr_per_cell->decrRef();
771 MEDCouplingFieldDiscretizationPerCell::MEDCouplingFieldDiscretizationPerCell(const MEDCouplingFieldDiscretizationPerCell& other, const int *startCellIds, const int *endCellIds):_discr_per_cell(0)
773 DataArrayInt *arr=other._discr_per_cell;
776 if(startCellIds==0 && endCellIds==0)
777 _discr_per_cell=arr->deepCpy();
779 _discr_per_cell=arr->selectByTupleIdSafe(startCellIds,endCellIds);
783 void MEDCouplingFieldDiscretizationPerCell::updateTime() const
786 updateTimeWith(*_discr_per_cell);
789 std::size_t MEDCouplingFieldDiscretizationPerCell::getHeapMemorySize() const
793 ret+=_discr_per_cell->getHeapMemorySize();
797 void MEDCouplingFieldDiscretizationPerCell::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception)
800 throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell has no discretization per cell !");
801 int nbOfTuples=_discr_per_cell->getNumberOfTuples();
802 if(nbOfTuples!=mesh->getNumberOfCells())
803 throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell has a discretization per cell but it's not matching the underlying mesh !");
806 bool MEDCouplingFieldDiscretizationPerCell::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
808 const MEDCouplingFieldDiscretizationPerCell *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationPerCell *>(other);
811 reason="Spatial discrtization of this is ON_GAUSS, which is not the case of other.";
814 if(_discr_per_cell==0)
815 return otherC->_discr_per_cell==0;
816 if(otherC->_discr_per_cell==0)
818 bool ret=_discr_per_cell->isEqualIfNotWhy(*otherC->_discr_per_cell,reason);
820 reason.insert(0,"Field discretization per cell DataArrayInt given the discid per cell :");
824 bool MEDCouplingFieldDiscretizationPerCell::isEqualWithoutConsideringStr(const MEDCouplingFieldDiscretization *other, double eps) const
826 const MEDCouplingFieldDiscretizationPerCell *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationPerCell *>(other);
829 if(_discr_per_cell==0)
830 return otherC->_discr_per_cell==0;
831 if(otherC->_discr_per_cell==0)
833 return _discr_per_cell->isEqualWithoutConsideringStr(*otherC->_discr_per_cell);
837 * This method is typically the first step of renumbering. The impact on _discr_per_cell is necessary here.
838 * virtualy by this method.
840 void MEDCouplingFieldDiscretizationPerCell::renumberCells(const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
842 int nbCells=_discr_per_cell->getNumberOfTuples();
843 const int *array=old2NewBg;
845 array=DataArrayInt::CheckAndPreparePermutation(old2NewBg,old2NewBg+nbCells);
847 DataArrayInt *dpc=_discr_per_cell->renumber(array);
848 _discr_per_cell->decrRef();
852 delete [] const_cast<int *>(array);
855 void MEDCouplingFieldDiscretizationPerCell::buildDiscrPerCellIfNecessary(const MEDCouplingMesh *m)
859 _discr_per_cell=DataArrayInt::New();
860 int nbTuples=m->getNumberOfCells();
861 _discr_per_cell->alloc(nbTuples,1);
862 int *ptr=_discr_per_cell->getPointer();
863 std::fill(ptr,ptr+nbTuples,DFT_INVALID_LOCID_VALUE);
867 void MEDCouplingFieldDiscretizationPerCell::checkNoOrphanCells() const throw(INTERP_KERNEL::Exception)
870 throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell::checkNoOrphanCells : no discretization defined !");
871 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> test=_discr_per_cell->getIdsEqual(DFT_INVALID_LOCID_VALUE);
872 if(test->getNumberOfTuples()!=0)
873 throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell::checkNoOrphanCells : presence of orphan cells !");
876 const DataArrayInt *MEDCouplingFieldDiscretizationPerCell::getArrayOfDiscIds() const
878 return _discr_per_cell;
881 MEDCouplingFieldDiscretizationGauss::MEDCouplingFieldDiscretizationGauss()
885 MEDCouplingFieldDiscretizationGauss::MEDCouplingFieldDiscretizationGauss(const MEDCouplingFieldDiscretizationGauss& other, const int *startCellIds, const int *endCellIds):MEDCouplingFieldDiscretizationPerCell(other,startCellIds,endCellIds),_loc(other._loc)
889 TypeOfField MEDCouplingFieldDiscretizationGauss::getEnum() const
894 bool MEDCouplingFieldDiscretizationGauss::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
896 const MEDCouplingFieldDiscretizationGauss *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationGauss *>(other);
899 reason="Spatial discrtization of this is ON_GAUSS, which is not the case of other.";
902 if(!MEDCouplingFieldDiscretizationPerCell::isEqualIfNotWhy(other,eps,reason))
904 if(_loc.size()!=otherC->_loc.size())
906 reason="Gauss spatial discretization : localization sizes differ";
909 std::size_t sz=_loc.size();
910 for(std::size_t i=0;i<sz;i++)
911 if(!_loc[i].isEqual(otherC->_loc[i],eps))
913 std::ostringstream oss; oss << "Gauss spatial discretization : Localization #" << i << " differ from this to other.";
920 bool MEDCouplingFieldDiscretizationGauss::isEqualWithoutConsideringStr(const MEDCouplingFieldDiscretization *other, double eps) const
922 const MEDCouplingFieldDiscretizationGauss *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationGauss *>(other);
925 if(!MEDCouplingFieldDiscretizationPerCell::isEqualWithoutConsideringStr(other,eps))
927 if(_loc.size()!=otherC->_loc.size())
929 std::size_t sz=_loc.size();
930 for(std::size_t i=0;i<sz;i++)
931 if(!_loc[i].isEqual(otherC->_loc[i],eps))
936 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationGauss::clone() const
938 return new MEDCouplingFieldDiscretizationGauss(*this);
941 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationGauss::clonePart(const int *startCellIds, const int *endCellIds) const
943 return new MEDCouplingFieldDiscretizationGauss(*this,startCellIds,endCellIds);
946 std::string MEDCouplingFieldDiscretizationGauss::getStringRepr() const
948 std::ostringstream oss; oss << REPR << "." << std::endl;
951 if(_discr_per_cell->isAllocated())
953 oss << "Discretization per cell : ";
954 std::copy(_discr_per_cell->begin(),_discr_per_cell->end(),std::ostream_iterator<int>(oss,", "));
958 oss << "Presence of " << _loc.size() << " localizations." << std::endl;
960 for(std::vector<MEDCouplingGaussLocalization>::const_iterator it=_loc.begin();it!=_loc.end();it++,i++)
962 oss << "+++++ Localization #" << i << " +++++" << std::endl;
963 oss << (*it).getStringRepr();
964 oss << "++++++++++" << std::endl;
969 std::size_t MEDCouplingFieldDiscretizationGauss::getHeapMemorySize() const
971 std::size_t ret=_loc.capacity()*sizeof(MEDCouplingGaussLocalization);
972 for(std::vector<MEDCouplingGaussLocalization>::const_iterator it=_loc.begin();it!=_loc.end();it++)
973 ret+=(*it).getHeapMemorySize();
974 return MEDCouplingFieldDiscretizationPerCell::getHeapMemorySize()+ret;
977 const char *MEDCouplingFieldDiscretizationGauss::getRepr() const
982 int MEDCouplingFieldDiscretizationGauss::getNumberOfTuples(const MEDCouplingMesh *) const
985 const int *dcPtr=_discr_per_cell->getConstPointer();
986 int nbOfTuples=_discr_per_cell->getNumberOfTuples();
987 for(const int *w=dcPtr;w!=dcPtr+nbOfTuples;w++)
988 ret+=_loc[*w].getNumberOfGaussPt();
992 int MEDCouplingFieldDiscretizationGauss::getNumberOfMeshPlaces(const MEDCouplingMesh *mesh) const
994 return mesh->getNumberOfCells();
997 DataArrayInt *MEDCouplingFieldDiscretizationGauss::getOffsetArr(const MEDCouplingMesh *mesh) const
999 int nbOfTuples=mesh->getNumberOfCells();
1000 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
1001 ret->alloc(nbOfTuples+1,1);
1002 int *retPtr=ret->getPointer();
1003 const int *start=_discr_per_cell->getConstPointer();
1004 if(_discr_per_cell->getNumberOfTuples()!=nbOfTuples)
1005 throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getOffsetArr : mismatch between the mesh and the discretization ids array length !");
1006 int maxPossible=(int)_loc.size();
1008 for(int i=0;i<nbOfTuples;i++,start++)
1010 if(*start>=0 && *start<maxPossible)
1011 retPtr[i+1]=retPtr[i]+_loc[*start].getNumberOfGaussPt();
1014 std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGauss::getOffsetArr : At position #" << i << " the locid = " << *start << " whereas it should be in [0," << maxPossible << ") !";
1015 throw INTERP_KERNEL::Exception(oss.str().c_str());
1021 void MEDCouplingFieldDiscretizationGauss::renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArrayDouble *>& arrays,
1022 const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
1024 const int *array=old2NewBg;
1026 array=DataArrayInt::CheckAndPreparePermutation(old2NewBg,old2NewBg+mesh->getNumberOfCells());
1027 int nbOfCells=_discr_per_cell->getNumberOfTuples();
1028 int nbOfTuples=getNumberOfTuples(0);
1029 const int *dcPtr=_discr_per_cell->getConstPointer();
1030 int *array2=new int[nbOfTuples];//stores the final conversion array old2New to give to arrays in renumberInPlace.
1031 int *array3=new int[nbOfCells];//store for each cell in present dcp array (already renumbered) the offset needed by each cell in new numbering.
1033 for(int i=1;i<nbOfCells;i++)
1034 array3[i]=array3[i-1]+_loc[dcPtr[i-1]].getNumberOfGaussPt();
1036 for(int i=0;i<nbOfCells;i++)
1038 int nbOfGaussPt=_loc[dcPtr[array[i]]].getNumberOfGaussPt();
1039 for(int k=0;k<nbOfGaussPt;k++,j++)
1040 array2[j]=array3[array[i]]+k;
1043 for(std::vector<DataArrayDouble *>::const_iterator it=arrays.begin();it!=arrays.end();it++)
1045 (*it)->renumberInPlace(array2);
1048 delete [] const_cast<int*>(array);
1051 DataArrayDouble *MEDCouplingFieldDiscretizationGauss::getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const
1053 checkNoOrphanCells();
1054 MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> umesh=mesh->buildUnstructured();//in general do nothing
1055 int nbOfTuples=getNumberOfTuples(mesh);
1056 DataArrayDouble *ret=DataArrayDouble::New();
1057 int spaceDim=mesh->getSpaceDimension();
1058 ret->alloc(nbOfTuples,spaceDim);
1059 std::vector< int > locIds;
1060 std::vector<DataArrayInt *> parts=splitIntoSingleGaussDicrPerCellType(locIds);
1061 std::vector< MEDCouplingAutoRefCountObjectPtr<DataArrayInt> > parts2(parts.size());
1062 std::copy(parts.begin(),parts.end(),parts2.begin());
1063 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> offsets=buildNbOfGaussPointPerCellField();
1064 offsets->computeOffsets();
1065 const int *ptrOffsets=offsets->getConstPointer();
1066 const double *coords=umesh->getCoords()->getConstPointer();
1067 const int *connI=umesh->getNodalConnectivityIndex()->getConstPointer();
1068 const int *conn=umesh->getNodalConnectivity()->getConstPointer();
1069 double *valsToFill=ret->getPointer();
1070 for(std::size_t i=0;i<parts2.size();i++)
1072 INTERP_KERNEL::GaussCoords calculator;
1074 const MEDCouplingGaussLocalization& cli=_loc[locIds[i]];//curLocInfo
1075 INTERP_KERNEL::NormalizedCellType typ=cli.getType();
1076 const std::vector<double>& wg=cli.getWeights();
1077 calculator.addGaussInfo(typ,INTERP_KERNEL::CellModel::GetCellModel(typ).getDimension(),
1078 &cli.getGaussCoords()[0],(int)wg.size(),&cli.getRefCoords()[0],
1079 INTERP_KERNEL::CellModel::GetCellModel(typ).getNumberOfNodes());
1081 int nbt=parts2[i]->getNumberOfTuples();
1082 for(const int *w=parts2[i]->getConstPointer();w!=parts2[i]->getConstPointer()+nbt;w++)
1083 calculator.calculateCoords(cli.getType(),coords,spaceDim,conn+connI[*w]+1,valsToFill+spaceDim*(ptrOffsets[*w]));
1085 ret->copyStringInfoFrom(*umesh->getCoords());
1089 void MEDCouplingFieldDiscretizationGauss::computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const int *partBg, const int *partEnd,
1090 DataArrayInt *&cellRest)
1092 throw INTERP_KERNEL::Exception("Not implemented yet !");
1098 void MEDCouplingFieldDiscretizationGauss::checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception)
1102 void MEDCouplingFieldDiscretizationGauss::getTinySerializationIntInformation(std::vector<int>& tinyInfo) const
1106 val=_discr_per_cell->getNumberOfTuples();
1107 tinyInfo.push_back(val);
1108 tinyInfo.push_back((int)_loc.size());
1110 tinyInfo.push_back(-1);
1112 tinyInfo.push_back(_loc[0].getDimension());
1113 for(std::vector<MEDCouplingGaussLocalization>::const_iterator iter=_loc.begin();iter!=_loc.end();iter++)
1114 (*iter).pushTinySerializationIntInfo(tinyInfo);
1117 void MEDCouplingFieldDiscretizationGauss::getTinySerializationDbleInformation(std::vector<double>& tinyInfo) const
1119 for(std::vector<MEDCouplingGaussLocalization>::const_iterator iter=_loc.begin();iter!=_loc.end();iter++)
1120 (*iter).pushTinySerializationDblInfo(tinyInfo);
1123 void MEDCouplingFieldDiscretizationGauss::getSerializationIntArray(DataArrayInt *& arr) const
1127 arr=_discr_per_cell;
1130 void MEDCouplingFieldDiscretizationGauss::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *& arr)
1132 int val=tinyInfo[0];
1135 _discr_per_cell=DataArrayInt::New();
1136 _discr_per_cell->alloc(val,1);
1140 arr=_discr_per_cell;
1141 int nbOfLoc=tinyInfo[1];
1143 int dim=tinyInfo[2];
1146 delta=((int)tinyInfo.size()-3)/nbOfLoc;
1147 for(int i=0;i<nbOfLoc;i++)
1149 std::vector<int> tmp(tinyInfo.begin()+3+i*delta,tinyInfo.begin()+3+(i+1)*delta);
1150 MEDCouplingGaussLocalization elt=MEDCouplingGaussLocalization::BuildNewInstanceFromTinyInfo(dim,tmp);
1151 _loc.push_back(elt);
1155 void MEDCouplingFieldDiscretizationGauss::finishUnserialization(const std::vector<double>& tinyInfo)
1157 double *tmp=new double[tinyInfo.size()];
1158 std::copy(tinyInfo.begin(),tinyInfo.end(),tmp);
1159 const double *work=tmp;
1160 for(std::vector<MEDCouplingGaussLocalization>::iterator iter=_loc.begin();iter!=_loc.end();iter++)
1161 work=(*iter).fillWithValues(work);
1165 double MEDCouplingFieldDiscretizationGauss::getIJK(const MEDCouplingMesh *mesh, const DataArrayDouble *da,
1166 int cellId, int nodeIdInCell, int compoId) const throw(INTERP_KERNEL::Exception)
1168 int offset=getOffsetOfCell(cellId);
1169 return da->getIJ(offset+nodeIdInCell,compoId);
1172 void MEDCouplingFieldDiscretizationGauss::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception)
1174 MEDCouplingFieldDiscretizationPerCell::checkCoherencyBetween(mesh,da);
1175 for(std::vector<MEDCouplingGaussLocalization>::const_iterator iter=_loc.begin();iter!=_loc.end();iter++)
1176 (*iter).checkCoherency();
1177 int nbOfDesc=(int)_loc.size();
1178 int nbOfCells=mesh->getNumberOfCells();
1179 const int *dc=_discr_per_cell->getConstPointer();
1180 for(int i=0;i<nbOfCells;i++)
1184 std::ostringstream oss; oss << "Cell # " << i << " of mesh \"" << mesh->getName() << "\" has an undefined gauss location ! Should never happend !";
1185 throw INTERP_KERNEL::Exception(oss.str().c_str());
1189 std::ostringstream oss; oss << "Cell # " << i << " of mesh \"" << mesh->getName() << "\" has no gauss location !";
1190 throw INTERP_KERNEL::Exception(oss.str().c_str());
1192 if(mesh->getTypeOfCell(i)!=_loc[dc[i]].getType())
1194 std::ostringstream oss; oss << "Types of mesh and gauss location mismatch for cell # " << i;
1195 throw INTERP_KERNEL::Exception(oss.str().c_str());
1198 int nbOfTuples=getNumberOfTuples(mesh);
1199 if(nbOfTuples!=da->getNumberOfTuples())
1201 std::ostringstream oss; oss << "Invalid number of tuples in the array : expecting " << nbOfTuples << " !";
1202 throw INTERP_KERNEL::Exception(oss.str().c_str());
1206 MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGauss::getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const
1209 throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getMeasureField : mesh instance specified is NULL !");
1210 throw INTERP_KERNEL::Exception("Not implemented yet !");
1213 void MEDCouplingFieldDiscretizationGauss::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
1215 throw INTERP_KERNEL::Exception("Not implemented yet !");
1218 void MEDCouplingFieldDiscretizationGauss::getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const
1220 throw INTERP_KERNEL::Exception("getValueOnPos(i,j,k) : Not applyable for Gauss points !");
1223 DataArrayDouble *MEDCouplingFieldDiscretizationGauss::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfPoints) const
1225 throw INTERP_KERNEL::Exception("getValueOnMulti : Not implemented yet for gauss points !");
1228 MEDCouplingMesh *MEDCouplingFieldDiscretizationGauss::buildSubMeshData(const MEDCouplingMesh *mesh, const int *start, const int *end, DataArrayInt *&di) const
1230 di=computeTupleIdsToSelectFromCellIds(mesh,start,end);
1231 return mesh->buildPart(start,end);
1235 * This method returns a tuple ids selection from cell ids selection [start;end).
1236 * This method is called by MEDCouplingFieldDiscretizationGauss::buildSubMeshData to return parameter \b di.
1238 * \return a newly allocated array containing ids to select into the DataArrayDouble of the field.
1241 DataArrayInt *MEDCouplingFieldDiscretizationGauss::computeTupleIdsToSelectFromCellIds(const MEDCouplingMesh *mesh, const int *startCellIds, const int *endCellIds) const
1244 throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::computeTupleIdsToSelectFromCellIds : null mesh !");
1245 if(!_discr_per_cell)
1246 throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::computeTupleIdsToSelectFromCellIds : null discretization ids !");
1247 int nbOfCells=mesh->getNumberOfCells();
1248 if(_discr_per_cell->getNumberOfTuples()!=nbOfCells)
1249 throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::computeTupleIdsToSelectFromCellIds : mismatch of nb of tuples of cell ids array and number of cells !");
1250 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfNodesPerCell=DataArrayInt::New(); nbOfNodesPerCell->alloc(nbOfCells,1);
1251 int *retPtr=nbOfNodesPerCell->getPointer();
1252 const int *pt=_discr_per_cell->getConstPointer();
1253 int nbMaxOfLocId=(int)_loc.size();
1254 for(int i=0;i<nbOfCells;i++,retPtr++,pt++)
1256 if(*pt>=0 && *pt<nbMaxOfLocId)
1257 *retPtr=_loc[*pt].getNumberOfGaussPt();
1259 nbOfNodesPerCell->computeOffsets2();
1260 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> sel=DataArrayInt::New(); sel->useArray(startCellIds,false,CPP_DEALLOC,(int)std::distance(startCellIds,endCellIds),1);
1261 return sel->buildExplicitArrByRanges(nbOfNodesPerCell);
1265 * No implementation needed !
1267 void MEDCouplingFieldDiscretizationGauss::renumberValuesOnNodes(double , const int *, int newNbOfNodes, DataArrayDouble *) const
1271 void MEDCouplingFieldDiscretizationGauss::renumberValuesOnCells(double epsOnVals, const MEDCouplingMesh *mesh, const int *old2New, int newSz, DataArrayDouble *arr) const
1273 throw INTERP_KERNEL::Exception("Not implemented yet !");
1276 void MEDCouplingFieldDiscretizationGauss::renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const
1278 throw INTERP_KERNEL::Exception("Number of cells has changed and becomes higher with some cells that have been split ! Unable to conserve the Gauss field !");
1281 void MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnType(const MEDCouplingMesh *m, INTERP_KERNEL::NormalizedCellType type, const std::vector<double>& refCoo,
1282 const std::vector<double>& gsCoo, const std::vector<double>& wg) throw(INTERP_KERNEL::Exception)
1284 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
1285 if((int)cm.getDimension()!=m->getMeshDimension())
1287 std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnType : mismatch of dimensions ! MeshDim==" << m->getMeshDimension();
1288 oss << " whereas Type '" << cm.getRepr() << "' has dimension " << cm.getDimension() << " !";
1289 throw INTERP_KERNEL::Exception(oss.str().c_str());
1291 buildDiscrPerCellIfNecessary(m);
1292 int id=(int)_loc.size();
1293 MEDCouplingGaussLocalization elt(type,refCoo,gsCoo,wg);
1294 _loc.push_back(elt);
1295 int *ptr=_discr_per_cell->getPointer();
1296 int nbCells=m->getNumberOfCells();
1297 for(int i=0;i<nbCells;i++)
1298 if(m->getTypeOfCell(i)==type)
1300 zipGaussLocalizations();
1303 void MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnCells(const MEDCouplingMesh *m, const int *begin, const int *end, const std::vector<double>& refCoo,
1304 const std::vector<double>& gsCoo, const std::vector<double>& wg) throw(INTERP_KERNEL::Exception)
1306 buildDiscrPerCellIfNecessary(m);
1307 if(std::distance(begin,end)<1)
1308 throw INTERP_KERNEL::Exception("Size of [begin,end) must be equal or greater than 1 !");
1309 INTERP_KERNEL::NormalizedCellType type=m->getTypeOfCell(*begin);
1310 MEDCouplingGaussLocalization elt(type,refCoo,gsCoo,wg);
1311 int id=(int)_loc.size();
1312 int *ptr=_discr_per_cell->getPointer();
1313 for(const int *w=begin+1;w!=end;w++)
1315 if(m->getTypeOfCell(*w)!=type)
1317 std::ostringstream oss; oss << "The cell with id " << *w << " has been detected to be incompatible in the [begin,end) array specified !";
1318 throw INTERP_KERNEL::Exception(oss.str().c_str());
1322 for(const int *w2=begin;w2!=end;w2++)
1325 _loc.push_back(elt);
1326 zipGaussLocalizations();
1329 void MEDCouplingFieldDiscretizationGauss::clearGaussLocalizations() throw(INTERP_KERNEL::Exception)
1333 _discr_per_cell->decrRef();
1339 MEDCouplingGaussLocalization& MEDCouplingFieldDiscretizationGauss::getGaussLocalization(int locId) throw(INTERP_KERNEL::Exception)
1341 checkLocalizationId(locId);
1345 int MEDCouplingFieldDiscretizationGauss::getNbOfGaussLocalization() const throw(INTERP_KERNEL::Exception)
1347 return (int)_loc.size();
1350 int MEDCouplingFieldDiscretizationGauss::getGaussLocalizationIdOfOneCell(int cellId) const throw(INTERP_KERNEL::Exception)
1352 if(!_discr_per_cell)
1353 throw INTERP_KERNEL::Exception("No Gauss localization still set !");
1354 int locId=_discr_per_cell->getConstPointer()[cellId];
1356 throw INTERP_KERNEL::Exception("No Gauss localization set for the specified cell !");
1360 int MEDCouplingFieldDiscretizationGauss::getGaussLocalizationIdOfOneType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception)
1362 std::set<int> ret=getGaussLocalizationIdsOfOneType(type);
1364 throw INTERP_KERNEL::Exception("No gauss discretization found for the specified type !");
1366 throw INTERP_KERNEL::Exception("Several gauss discretizations have been found for the specified type !");
1367 return *ret.begin();
1370 std::set<int> MEDCouplingFieldDiscretizationGauss::getGaussLocalizationIdsOfOneType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception)
1372 if(!_discr_per_cell)
1373 throw INTERP_KERNEL::Exception("No Gauss localization still set !");
1376 for(std::vector<MEDCouplingGaussLocalization>::const_iterator iter=_loc.begin();iter!=_loc.end();iter++,id++)
1377 if((*iter).getType()==type)
1382 void MEDCouplingFieldDiscretizationGauss::getCellIdsHavingGaussLocalization(int locId, std::vector<int>& cellIds) const throw(INTERP_KERNEL::Exception)
1384 if(locId<0 || locId>=(int)_loc.size())
1385 throw INTERP_KERNEL::Exception("Invalid locId given : must be in range [0:getNbOfGaussLocalization()) !");
1386 int nbOfTuples=_discr_per_cell->getNumberOfTuples();
1387 const int *ptr=_discr_per_cell->getConstPointer();
1388 for(int i=0;i<nbOfTuples;i++)
1390 cellIds.push_back(i);
1393 const MEDCouplingGaussLocalization& MEDCouplingFieldDiscretizationGauss::getGaussLocalization(int locId) const throw(INTERP_KERNEL::Exception)
1395 checkLocalizationId(locId);
1399 void MEDCouplingFieldDiscretizationGauss::checkLocalizationId(int locId) const throw(INTERP_KERNEL::Exception)
1401 if(locId<0 || locId>=(int)_loc.size())
1402 throw INTERP_KERNEL::Exception("Invalid locId given : must be in range [0:getNbOfGaussLocalization()) !");
1405 int MEDCouplingFieldDiscretizationGauss::getOffsetOfCell(int cellId) const throw(INTERP_KERNEL::Exception)
1408 const int *start=_discr_per_cell->getConstPointer();
1409 for(const int *w=start;w!=start+cellId;w++)
1410 ret+=_loc[*w].getNumberOfGaussPt();
1415 * This method do the assumption that there is no orphan cell. If there is an exception is thrown.
1416 * This method makes the assumption too that '_discr_per_cell' is defined. If not an exception is thrown.
1417 * This method returns a newly created array with number of tuples equals to '_discr_per_cell->getNumberOfTuples' and number of components equal to 1.
1418 * The i_th tuple in returned array is the number of gauss point if the corresponding cell.
1420 DataArrayInt *MEDCouplingFieldDiscretizationGauss::buildNbOfGaussPointPerCellField() const throw(INTERP_KERNEL::Exception)
1422 if(!_discr_per_cell)
1423 throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::buildNbOfGaussPointPerCellField : no discretization array set !");
1424 int nbOfTuples=_discr_per_cell->getNumberOfTuples();
1425 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
1426 const int *w=_discr_per_cell->getConstPointer();
1427 ret->alloc(nbOfTuples,1);
1428 int *valsToFill=ret->getPointer();
1429 for(int i=0;i<nbOfTuples;i++,w++)
1430 if(*w!=DFT_INVALID_LOCID_VALUE)
1431 valsToFill[i]=_loc[*w].getNumberOfGaussPt();
1433 throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::buildNbOfGaussPointPerCellField : orphan cell detected !");
1438 * This method makes the assumption that _discr_per_cell is set.
1439 * This method reduces as much as possible number size of _loc.
1440 * This method is usefull when several set on same cells has been done and that some Gauss Localization are no more used.
1442 void MEDCouplingFieldDiscretizationGauss::zipGaussLocalizations()
1444 const int *start=_discr_per_cell->getConstPointer();
1445 int nbOfTuples=_discr_per_cell->getNumberOfTuples();
1446 INTERP_KERNEL::AutoPtr<int> tmp=new int[_loc.size()];
1447 std::fill((int *)tmp,(int *)tmp+_loc.size(),-2);
1448 for(const int *w=start;w!=start+nbOfTuples;w++)
1452 for(int i=0;i<(int)_loc.size();i++)
1455 if(fid==(int)_loc.size())
1458 int *start2=_discr_per_cell->getPointer();
1459 for(int *w2=start2;w2!=start2+nbOfTuples;w2++)
1462 std::vector<MEDCouplingGaussLocalization> tmpLoc;
1463 for(int i=0;i<(int)_loc.size();i++)
1465 tmpLoc.push_back(_loc[tmp[i]]);
1470 * This method is usefull when 'this' describes a field discretization with several gauss discretization on a \b same cell type.
1471 * For example same NORM_TRI3 cells having 6 gauss points and others with 12 gauss points.
1472 * This method returns 2 arrays with same size : the return value and 'locIds' output parameter.
1473 * For a given i into [0,locIds.size) ret[i] represents the set of cell ids of i_th set an locIds[i] represents the set of discretisation of the set.
1474 * The return vector contains a set of newly created instance to deal with.
1475 * The returned vector represents a \b partition of cells ids with a gauss discretization set.
1477 * If no descretization is set in 'this' and exception will be thrown.
1479 std::vector<DataArrayInt *> MEDCouplingFieldDiscretizationGauss::splitIntoSingleGaussDicrPerCellType(std::vector<int>& locIds) const throw(INTERP_KERNEL::Exception)
1481 if(!_discr_per_cell)
1482 throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::splitIntoSingleGaussDicrPerCellType : no descretization set !");
1483 return _discr_per_cell->partitionByDifferentValues(locIds);
1486 MEDCouplingFieldDiscretizationGaussNE::MEDCouplingFieldDiscretizationGaussNE()
1490 TypeOfField MEDCouplingFieldDiscretizationGaussNE::getEnum() const
1495 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationGaussNE::clone() const
1497 return new MEDCouplingFieldDiscretizationGaussNE(*this);
1500 std::string MEDCouplingFieldDiscretizationGaussNE::getStringRepr() const
1502 return std::string(REPR);
1505 const char *MEDCouplingFieldDiscretizationGaussNE::getRepr() const
1510 bool MEDCouplingFieldDiscretizationGaussNE::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
1512 const MEDCouplingFieldDiscretizationGaussNE *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationGaussNE *>(other);
1515 reason="Spatial discrtization of this is ON_GAUSS_NE, which is not the case of other.";
1519 int MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuples(const MEDCouplingMesh *mesh) const
1522 int nbOfCells=mesh->getNumberOfCells();
1523 for(int i=0;i<nbOfCells;i++)
1525 INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell(i);
1526 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
1528 throw INTERP_KERNEL::Exception("Not implemented yet Gauss node on elements for polygons and polyedrons !");
1529 ret+=cm.getNumberOfNodes();
1534 int MEDCouplingFieldDiscretizationGaussNE::getNumberOfMeshPlaces(const MEDCouplingMesh *mesh) const
1536 return mesh->getNumberOfCells();
1539 DataArrayInt *MEDCouplingFieldDiscretizationGaussNE::getOffsetArr(const MEDCouplingMesh *mesh) const
1541 int nbOfTuples=mesh->getNumberOfCells();
1542 DataArrayInt *ret=DataArrayInt::New();
1543 ret->alloc(nbOfTuples+1,1);
1544 int *retPtr=ret->getPointer();
1546 for(int i=0;i<nbOfTuples;i++)
1548 INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell(i);
1549 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
1551 throw INTERP_KERNEL::Exception("Not implemented yet Gauss node on elements for polygons and polyedrons !");
1552 retPtr[i+1]=retPtr[i]+cm.getNumberOfNodes();
1557 void MEDCouplingFieldDiscretizationGaussNE::renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArrayDouble *>& arrays,
1558 const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
1560 const int *array=old2NewBg;
1562 array=DataArrayInt::CheckAndPreparePermutation(old2NewBg,old2NewBg+mesh->getNumberOfCells());
1563 int nbOfCells=mesh->getNumberOfCells();
1564 int nbOfTuples=getNumberOfTuples(mesh);
1565 int *array2=new int[nbOfTuples];//stores the final conversion array old2New to give to arrays in renumberInPlace.
1566 int *array3=new int[nbOfCells];//store for each cell in after renumbering the offset needed by each cell in new numbering.
1568 for(int i=1;i<nbOfCells;i++)
1570 INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell((int)std::distance(array,std::find(array,array+nbOfCells,i-1)));
1571 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
1572 array3[i]=array3[i-1]+cm.getNumberOfNodes();
1575 for(int i=0;i<nbOfCells;i++)
1577 INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell(i);
1578 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
1579 for(int k=0;k<(int)cm.getNumberOfNodes();k++,j++)
1580 array2[j]=array3[array[i]]+k;
1583 for(std::vector<DataArrayDouble *>::const_iterator it=arrays.begin();it!=arrays.end();it++)
1585 (*it)->renumberInPlace(array2);
1588 delete [] const_cast<int *>(array);
1591 DataArrayDouble *MEDCouplingFieldDiscretizationGaussNE::getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const
1593 throw INTERP_KERNEL::Exception("Not implemented yet !");
1596 void MEDCouplingFieldDiscretizationGaussNE::computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const int *partBg, const int *partEnd,
1597 DataArrayInt *&cellRest)
1599 throw INTERP_KERNEL::Exception("Not implemented yet !");
1602 void MEDCouplingFieldDiscretizationGaussNE::checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception)
1606 double MEDCouplingFieldDiscretizationGaussNE::getIJK(const MEDCouplingMesh *mesh, const DataArrayDouble *da,
1607 int cellId, int nodeIdInCell, int compoId) const throw(INTERP_KERNEL::Exception)
1610 for(int i=0;i<cellId;i++)
1612 INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell(i);
1613 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
1614 offset+=cm.getNumberOfNodes();
1616 return da->getIJ(offset+nodeIdInCell,compoId);
1619 void MEDCouplingFieldDiscretizationGaussNE::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception)
1621 int nbOfTuples=getNumberOfTuples(mesh);
1622 if(nbOfTuples!=da->getNumberOfTuples())
1624 std::ostringstream oss; oss << "Invalid number of tuples in the array : expecting " << nbOfTuples << " !";
1625 throw INTERP_KERNEL::Exception(oss.str().c_str());
1629 MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGaussNE::getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const
1632 throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::getMeasureField : mesh instance specified is NULL !");
1633 throw INTERP_KERNEL::Exception("Not implemented yet !");
1636 void MEDCouplingFieldDiscretizationGaussNE::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
1638 throw INTERP_KERNEL::Exception("Not implemented yet !");
1641 void MEDCouplingFieldDiscretizationGaussNE::getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const
1643 throw INTERP_KERNEL::Exception("getValueOnPos(i,j,k) : Not applyable for Gauss points !");
1646 DataArrayDouble *MEDCouplingFieldDiscretizationGaussNE::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfPoints) const
1648 throw INTERP_KERNEL::Exception("getValueOnMulti : Not implemented for Gauss NE !");
1651 MEDCouplingMesh *MEDCouplingFieldDiscretizationGaussNE::buildSubMeshData(const MEDCouplingMesh *mesh, const int *start, const int *end, DataArrayInt *&di) const
1653 di=computeTupleIdsToSelectFromCellIds(mesh,start,end);
1654 return mesh->buildPart(start,end);
1658 * This method returns a tuple ids selection from cell ids selection [start;end).
1659 * This method is called by MEDCouplingFieldDiscretizationGaussNE::buildSubMeshData to return parameter \b di.
1661 * \return a newly allocated array containing ids to select into the DataArrayDouble of the field.
1664 DataArrayInt *MEDCouplingFieldDiscretizationGaussNE::computeTupleIdsToSelectFromCellIds(const MEDCouplingMesh *mesh, const int *startCellIds, const int *endCellIds) const
1667 throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::computeTupleIdsToSelectFromCellIds : null mesh !");
1668 const MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> umesh=mesh->buildUnstructured();
1669 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfNodesPerCell=umesh->computeNbOfNodesPerCell();
1670 nbOfNodesPerCell->computeOffsets2();
1671 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> sel=DataArrayInt::New(); sel->useArray(startCellIds,false,CPP_DEALLOC,(int)std::distance(startCellIds,endCellIds),1);
1672 return sel->buildExplicitArrByRanges(nbOfNodesPerCell);
1676 * No implementation needed !
1678 void MEDCouplingFieldDiscretizationGaussNE::renumberValuesOnNodes(double , const int *, int newNbOfNodes, DataArrayDouble *) const
1682 void MEDCouplingFieldDiscretizationGaussNE::renumberValuesOnCells(double epsOnVals, const MEDCouplingMesh *mesh, const int *old2New, int newSz, DataArrayDouble *arr) const
1684 throw INTERP_KERNEL::Exception("Not implemented yet !");
1687 void MEDCouplingFieldDiscretizationGaussNE::renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const
1689 throw INTERP_KERNEL::Exception("Not implemented yet !");
1692 MEDCouplingFieldDiscretizationGaussNE::MEDCouplingFieldDiscretizationGaussNE(const MEDCouplingFieldDiscretizationGaussNE& other):MEDCouplingFieldDiscretization(other)
1696 TypeOfField MEDCouplingFieldDiscretizationKriging::getEnum() const
1701 const char *MEDCouplingFieldDiscretizationKriging::getRepr() const
1706 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationKriging::clone() const
1708 return new MEDCouplingFieldDiscretizationKriging;
1711 std::string MEDCouplingFieldDiscretizationKriging::getStringRepr() const
1713 return std::string(REPR);
1716 void MEDCouplingFieldDiscretizationKriging::checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception)
1718 if(nat!=ConservativeVolumic)
1719 throw INTERP_KERNEL::Exception("Invalid nature for Kriging field : expected ConservativeVolumic !");
1722 bool MEDCouplingFieldDiscretizationKriging::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
1724 const MEDCouplingFieldDiscretizationKriging *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationKriging *>(other);
1727 reason="Spatial discrtization of this is ON_NODES_KR, which is not the case of other.";
1731 MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationKriging::getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const
1734 throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::getMeasureField : mesh instance specified is NULL !");
1735 throw INTERP_KERNEL::Exception("getMeasureField on FieldDiscretizationKriging : not implemented yet !");
1738 void MEDCouplingFieldDiscretizationKriging::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
1740 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> res2=MEDCouplingFieldDiscretizationKriging::getValueOnMulti(arr,mesh,loc,1);
1741 std::copy(res2->begin(),res2->end(),res);
1744 DataArrayDouble *MEDCouplingFieldDiscretizationKriging::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfTargetPoints) const
1746 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords=getLocalizationOfDiscValues(mesh);
1747 int nbOfPts=coords->getNumberOfTuples();
1748 int dimension=coords->getNumberOfComponents();
1751 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> KnewiK=computeVectorOfCoefficients(mesh,arr,delta);
1753 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> locArr=DataArrayDouble::New();
1754 locArr->useArray(loc,false,CPP_DEALLOC,nbOfTargetPoints,dimension);
1755 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrix2=coords->buildEuclidianDistanceDenseMatrixWith(locArr);
1756 operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfPts*nbOfTargetPoints,matrix2->getPointer());
1757 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrix3=DataArrayDouble::New();
1758 matrix3->alloc((nbOfPts+delta)*nbOfTargetPoints,1);
1759 double *work=matrix3->getPointer();
1760 const double *workCst=matrix2->getConstPointer();
1761 const double *workCst2=loc;
1762 for(int i=0;i<nbOfTargetPoints;i++,workCst+=nbOfPts,workCst2+=delta-1)
1764 for(int j=0;j<nbOfPts;j++)
1765 work[j*nbOfTargetPoints+i]=workCst[j];
1766 work[nbOfPts*nbOfTargetPoints+i]=1.0;
1767 for(int j=0;j<delta-1;j++)
1768 work[(nbOfPts+1+j)*nbOfTargetPoints+i]=workCst2[j];
1771 int nbOfCompo=arr->getNumberOfComponents();
1772 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=DataArrayDouble::New();
1773 ret->alloc(nbOfTargetPoints,nbOfCompo);
1774 INTERP_KERNEL::matrixProduct(KnewiK->getConstPointer(),1,nbOfPts+delta,matrix3->getConstPointer(),nbOfPts+delta,nbOfTargetPoints*nbOfCompo,ret->getPointer());
1779 * This method computes coefficients to apply to each representing points of \a mesh, that is to say the nodes of \a mesh given a field array \a arr whose
1780 * number of tuples should be equal to the number of representing points in \a mesh.
1782 * \param [in] mesh is the sources of nodes on which kriging will be done regarding the parameters and the value of \c this->getSpaceDimension()
1783 * \param [in] arr input field DataArrayDouble whose number of tuples must be equal to the number of nodes in \a mesh
1784 * \param [out] isDrift return if drift coefficients are present in the returned vector of coefficients, and if. If different from 0 there is presence of drift coefficients.
1785 * Whatever the value of \a isDrift the number of tuples of returned DataArrayDouble will be equal to \c arr->getNumberOfTuples() + \a isDrift.
1786 * \return a newly allocated array containing coefficients including or not drift coefficient at the end depending the value of \a isDrift parameter.
1788 DataArrayDouble *MEDCouplingFieldDiscretizationKriging::computeVectorOfCoefficients(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, int& isDrift) const
1790 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords=getLocalizationOfDiscValues(mesh);
1791 int nbOfPts=coords->getNumberOfTuples();
1792 //int dimension=coords->getNumberOfComponents();
1793 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrix=coords->buildEuclidianDistanceDenseMatrix();
1794 operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfPts*nbOfPts,matrix->getPointer());
1796 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixWithDrift=performDrift(matrix,coords,isDrift);
1797 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixInv=DataArrayDouble::New();
1798 matrixInv->alloc((nbOfPts+isDrift)*(nbOfPts+isDrift),1);
1799 INTERP_KERNEL::inverseMatrix(matrixWithDrift->getConstPointer(),nbOfPts+isDrift,matrixInv->getPointer());
1801 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> KnewiK=DataArrayDouble::New();
1802 KnewiK->alloc((nbOfPts+isDrift)*1,1);
1803 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr2=DataArrayDouble::New();
1804 arr2->alloc((nbOfPts+isDrift)*1,1);
1805 double *work=std::copy(arr->begin(),arr->end(),arr2->getPointer());
1806 std::fill(work,work+isDrift,0.);
1807 INTERP_KERNEL::matrixProduct(matrixInv->getConstPointer(),nbOfPts+isDrift,nbOfPts+isDrift,arr2->getConstPointer(),nbOfPts+isDrift,1,KnewiK->getPointer());
1808 return KnewiK.retn();
1812 * Apply \f f(x) on each element x in \a matrixPtr. \a matrixPtr is expected to be a dense matrix represented by a chunck of memory of size at least equal to \a nbOfElems.
1814 * \param [in] spaceDimension space dimension of the input mesh on which the Kriging has to be performed
1815 * \param [in] nbOfElems is the result of the product of nb of rows and the nb of columns of matrix \a matrixPtr
1816 * \param [in,out] matrixPtr is the dense matrix whose on each values the operation will be applied
1818 void MEDCouplingFieldDiscretizationKriging::operateOnDenseMatrix(int spaceDimension, int nbOfElems, double *matrixPtr) const
1820 switch(spaceDimension)
1824 for(int i=0;i<nbOfElems;i++)
1826 double val=matrixPtr[i];
1827 matrixPtr[i]=val*val*val;
1832 throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::operateOnDenseMatrix : only dimension 1 implemented !");
1837 * Starting from a square matrix \a matr, this method returns a newly allocated dense square matrix whose \a matr is included in returned matrix
1838 * in the top left corner, and in the remaining returned matrix the parameters to take into account about the kriging drift.
1839 * For the moment only linear srift is implemented.
1841 * \param [in] arr the position of points were input mesh geometry is considered for Kriging
1842 * \param [in] matr input matrix whose drift part will be added
1843 * \param [out] delta the difference between the size of the output matrix and the input matrix \a matr.
1844 * \return a newly allocated matrix bigger than input matrix \a matr.
1846 DataArrayDouble *MEDCouplingFieldDiscretizationKriging::performDrift(const DataArrayDouble *matr, const DataArrayDouble *arr, int& delta) const
1848 int spaceDimension=arr->getNumberOfComponents();
1849 delta=spaceDimension+1;
1850 int szOfMatrix=arr->getNumberOfTuples();
1851 if(szOfMatrix*szOfMatrix!=matr->getNumberOfTuples())
1852 throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::performDrift : invalid size");
1853 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=DataArrayDouble::New();
1854 ret->alloc((szOfMatrix+delta)*(szOfMatrix+delta),1);
1855 const double *srcWork=matr->getConstPointer();
1856 const double *srcWork2=arr->getConstPointer();
1857 double *destWork=ret->getPointer();
1858 for(int i=0;i<szOfMatrix;i++)
1860 destWork=std::copy(srcWork,srcWork+szOfMatrix,destWork);
1861 srcWork+=szOfMatrix;
1863 destWork=std::copy(srcWork2,srcWork2+spaceDimension,destWork);
1864 srcWork2+=spaceDimension;
1866 std::fill(destWork,destWork+szOfMatrix,1.); destWork+=szOfMatrix;
1867 std::fill(destWork,destWork+spaceDimension+1,0.); destWork+=spaceDimension+1;
1868 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arrNoI=arr->toNoInterlace();
1869 srcWork2=arrNoI->getConstPointer();
1870 for(int i=0;i<spaceDimension;i++)
1872 destWork=std::copy(srcWork2,srcWork2+szOfMatrix,destWork);
1873 srcWork2+=szOfMatrix;
1874 std::fill(destWork,destWork+spaceDimension+1,0.);
1875 destWork+=spaceDimension;