1 // Copyright (C) 2007-2016 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, or (at your option) any later version.
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 "MEDCouplingGaussLocalization.hxx"
22 #include "InterpKernelGaussCoords.hxx"
23 #include "MEDCouplingUMesh.hxx"
24 #include "CellModel.hxx"
32 using namespace MEDCoupling;
34 MEDCouplingGaussLocalization::MEDCouplingGaussLocalization(INTERP_KERNEL::NormalizedCellType type, const std::vector<double>& refCoo,
35 const std::vector<double>& gsCoo, const std::vector<double>& w)
36 try:_type(type),_ref_coord(refCoo),_gauss_coord(gsCoo),_weight(w)
38 checkConsistencyLight();
40 catch(INTERP_KERNEL::Exception& e)
42 _type=INTERP_KERNEL::NORM_ERROR;
49 MEDCouplingGaussLocalization::MEDCouplingGaussLocalization(INTERP_KERNEL::NormalizedCellType typ)
52 INTERP_KERNEL::CellModel::GetCellModel(_type);
54 catch(INTERP_KERNEL::Exception& e)
56 _type=INTERP_KERNEL::NORM_ERROR;
60 void MEDCouplingGaussLocalization::setType(INTERP_KERNEL::NormalizedCellType typ)
62 INTERP_KERNEL::CellModel::GetCellModel(typ);//throws if not found. This is a check
66 void MEDCouplingGaussLocalization::checkConsistencyLight() const
68 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(_type);
69 int nbNodes=cm.getNumberOfNodes();
70 int dim=cm.getDimension();
73 if((int)_ref_coord.size()!=nbNodes*dim)
75 std::ostringstream oss; oss << "Invalid size of refCoo : expecting to be : " << nbNodes << " (nbNodePerCell) * " << dim << " (dim) !";
76 throw INTERP_KERNEL::Exception(oss.str().c_str());
79 if(_gauss_coord.size()!=dim*_weight.size())
81 std::ostringstream oss; oss << "Invalid gsCoo size and weight size : gsCoo.size() must be equal to _weight.size() * " << dim << " (dim) !";
82 throw INTERP_KERNEL::Exception(oss.str().c_str());
86 int MEDCouplingGaussLocalization::getDimension() const
90 return (int)_gauss_coord.size()/(int)_weight.size();
93 int MEDCouplingGaussLocalization::getNumberOfPtsInRefCell() const
95 int dim=getDimension();
98 return (int)_ref_coord.size()/dim;
101 std::string MEDCouplingGaussLocalization::getStringRepr() const
103 std::ostringstream oss;
104 oss << "CellType : " << INTERP_KERNEL::CellModel::GetCellModel(_type).getRepr() << std::endl;
105 oss << "Ref coords : "; std::copy(_ref_coord.begin(),_ref_coord.end(),std::ostream_iterator<double>(oss,", ")); oss << std::endl;
106 oss << "Localization coords : "; std::copy(_gauss_coord.begin(),_gauss_coord.end(),std::ostream_iterator<double>(oss,", ")); oss << std::endl;
107 oss << "Weight : "; std::copy(_weight.begin(),_weight.end(),std::ostream_iterator<double>(oss,", ")); oss << std::endl;
111 std::size_t MEDCouplingGaussLocalization::getMemorySize() const
114 ret+=_ref_coord.capacity()*sizeof(double);
115 ret+=_gauss_coord.capacity()*sizeof(double);
116 ret+=_weight.capacity()*sizeof(double);
120 bool MEDCouplingGaussLocalization::isEqual(const MEDCouplingGaussLocalization& other, double eps) const
122 if(_type!=other._type)
124 if(!AreAlmostEqual(_ref_coord,other._ref_coord,eps))
126 if(!AreAlmostEqual(_gauss_coord,other._gauss_coord,eps))
128 if(!AreAlmostEqual(_weight,other._weight,eps))
133 double MEDCouplingGaussLocalization::getRefCoord(int ptIdInCell, int comp) const
135 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(_type);
136 int nbNodes=cm.getNumberOfNodes();
137 int dim=cm.getDimension();
138 if(ptIdInCell<0 || ptIdInCell>=nbNodes)
139 throw INTERP_KERNEL::Exception("ptIdInCell specified is invalid : must be in [0;nbNodesPerCell) !");
140 if(comp<0 || comp>=dim)
141 throw INTERP_KERNEL::Exception("comp specified is invalid : must be in [0:dimOfCell) !");
142 return _ref_coord[ptIdInCell*dim+comp];
145 double MEDCouplingGaussLocalization::getGaussCoord(int gaussPtIdInCell, int comp) const
147 int dim=checkCoherencyOfRequest(gaussPtIdInCell,comp);
148 return _gauss_coord[gaussPtIdInCell*dim+comp];
151 double MEDCouplingGaussLocalization::getWeight(int gaussPtIdInCell, double newVal) const
153 checkCoherencyOfRequest(gaussPtIdInCell,0);
154 return _weight[gaussPtIdInCell];
158 * Completely useless method for end user. Only for CORBA MPI serialization/unserialization.
159 * push at the end of tinyInfo its basic serialization info. The size of pushed data is always the same.
160 * @param tinyInfo inout parameter.
162 void MEDCouplingGaussLocalization::pushTinySerializationIntInfo(std::vector<int>& tinyInfo) const
164 tinyInfo.push_back((int)_type);
165 tinyInfo.push_back(getNumberOfPtsInRefCell());
166 tinyInfo.push_back(getNumberOfGaussPt());
170 * Completely useless method for end user. Only for CORBA MPI serialization/unserialization.
171 * push at the end of tinyInfo its basic serialization info. The size of pushed data is \b NOT always the same contrary to pushTinySerializationIntInfo.
172 * @param tinyInfo inout parameter.
174 void MEDCouplingGaussLocalization::pushTinySerializationDblInfo(std::vector<double>& tinyInfo) const
176 tinyInfo.insert(tinyInfo.end(),_ref_coord.begin(),_ref_coord.end());
177 tinyInfo.insert(tinyInfo.end(),_gauss_coord.begin(),_gauss_coord.end());
178 tinyInfo.insert(tinyInfo.end(),_weight.begin(),_weight.end());
182 * This method operates the exact inverse operation than MEDCouplingGaussLocalization::pushTinySerializationDblInfo method. This is one of the last step of unserialization process.
183 * This method should be called on an object resized by buildNewInstanceFromTinyInfo static method.
184 * This method takes in argument a pointer 'vals' that point to the begin of double data pushed remotely by pushTinySerializationDblInfo method.
185 * This method returns the pointer 'vals' with an offset of size what it has been read in this method.
187 const double *MEDCouplingGaussLocalization::fillWithValues(const double *vals)
189 const double *work=vals;
190 std::copy(work,work+_ref_coord.size(),_ref_coord.begin());
191 work+=_ref_coord.size();
192 std::copy(work,work+_gauss_coord.size(),_gauss_coord.begin());
193 work+=_gauss_coord.size();
194 std::copy(work,work+_weight.size(),_weight.begin());
195 work+=_weight.size();
200 * Given points in \a ptsInRefCoo in the reference coordinates of \a this (in _ref_coord attribute)
201 * this method computes their coordinates in real world for each cells in \a mesh.
202 * So the returned array will contain nbCells* \a ptsInRefCoo->getNumberOfTuples() tuples and the number of component
203 * will be equal to the dimension of \a this.
205 * This method ignores Gauss points in \a this and only those in \a ptsInRefCoo are considered here.
207 MCAuto<DataArrayDouble> MEDCouplingGaussLocalization::localizePtsInRefCooForEachCell(const DataArrayDouble *ptsInRefCoo, const MEDCouplingUMesh *mesh) const
209 if(!ptsInRefCoo || !mesh)
210 throw INTERP_KERNEL::Exception("MEDCouplingGaussLocalization::localizePtsInRefCooForEachCell : null pointer !");
211 ptsInRefCoo->checkAllocated();
212 mesh->checkConsistencyLight();
214 int nbCells(mesh->getNumberOfCells());
215 const double *coords(mesh->getCoords()->begin());
216 const int *connI(mesh->getNodalConnectivityIndex()->begin()),*conn(mesh->getNodalConnectivity()->begin());
218 int nbPts(ptsInRefCoo->getNumberOfTuples());
219 INTERP_KERNEL::NormalizedCellType typ(getType());
220 int dim(INTERP_KERNEL::CellModel::GetCellModel(typ).getDimension());
221 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
222 ret->alloc(nbPts*nbCells,dim);
223 double *retPtr(ret->getPointer());
224 if(dim!=ptsInRefCoo->getNumberOfComponents())
225 throw INTERP_KERNEL::Exception("MEDCouplingGaussLocalization::localizePtsInRefCooForEachCell : number of components of input coo is not equal to dim of element !");
226 const std::vector<double>& wg(getWeights());
227 INTERP_KERNEL::GaussCoords calculator;
228 calculator.addGaussInfo(typ,dim, ptsInRefCoo->begin(),nbPts,&_ref_coord[0],getNumberOfPtsInRefCell());
230 for(int i=0;i<nbCells;i++,retPtr+=nbPts*dim)
231 calculator.calculateCoords(getType(),coords,dim,conn+connI[i]+1,retPtr);
236 * This method sets the comp_th component of ptIdInCell_th point coordinate of reference element of type this->_type.
237 * @throw if not 0<=ptIdInCell<nbOfNodePerCell or if not 0<=comp<dim
239 void MEDCouplingGaussLocalization::setRefCoord(int ptIdInCell, int comp, double newVal)
241 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(_type);
242 int nbNodes=cm.getNumberOfNodes();
243 int dim=cm.getDimension();
244 if(ptIdInCell<0 || ptIdInCell>=nbNodes)
245 throw INTERP_KERNEL::Exception("ptIdInCell specified is invalid : must be in [0;nbNodesPerCell) !");
246 if(comp<0 || comp>=dim)
247 throw INTERP_KERNEL::Exception("comp specified is invalid : must be in [0:dimOfCell) !");
248 _ref_coord[ptIdInCell*dim+comp]=newVal;
251 void MEDCouplingGaussLocalization::setGaussCoord(int gaussPtIdInCell, int comp, double newVal)
253 int dim=checkCoherencyOfRequest(gaussPtIdInCell,comp);
254 _gauss_coord[gaussPtIdInCell*dim+comp]=newVal;
257 void MEDCouplingGaussLocalization::setWeight(int gaussPtIdInCell, double newVal)
259 checkCoherencyOfRequest(gaussPtIdInCell,0);
260 _weight[gaussPtIdInCell]=newVal;
263 void MEDCouplingGaussLocalization::setRefCoords(const std::vector<double>& refCoo)
268 void MEDCouplingGaussLocalization::setGaussCoords(const std::vector<double>& gsCoo)
273 void MEDCouplingGaussLocalization::setWeights(const std::vector<double>& w)
279 * The format of 'tinyData' parameter is the same than pushed in method MEDCouplingGaussLocalization::pushTinySerializationIntInfo.
281 MEDCouplingGaussLocalization MEDCouplingGaussLocalization::BuildNewInstanceFromTinyInfo(int dim, const std::vector<int>& tinyData)
283 std::vector<double> v1(dim*tinyData[1]),v2(dim*tinyData[2]),v3(tinyData[2]);
284 return MEDCouplingGaussLocalization((INTERP_KERNEL::NormalizedCellType)tinyData[0],v1,v2,v3);
287 int MEDCouplingGaussLocalization::checkCoherencyOfRequest(int gaussPtIdInCell, int comp) const
289 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(_type);
290 int dim=cm.getDimension();
291 int nbGsPts=getNumberOfGaussPt();
292 if(gaussPtIdInCell<0 || gaussPtIdInCell>=nbGsPts)
293 throw INTERP_KERNEL::Exception("gaussPtIdInCell specified is invalid : must be in [0:nbGsPts) !");
294 if(comp<0 || comp>=dim)
295 throw INTERP_KERNEL::Exception("comp specified is invalid : must be in [0:dimOfCell) !");
299 bool MEDCouplingGaussLocalization::AreAlmostEqual(const std::vector<double>& v1, const std::vector<double>& v2, double eps)
301 std::size_t sz=v1.size();
304 std::vector<double> tmp(sz);
305 std::transform(v1.begin(),v1.end(),v2.begin(),tmp.begin(),std::minus<double>());
306 std::transform(tmp.begin(),tmp.end(),tmp.begin(),std::ptr_fun<double,double>(fabs));
307 return *std::max_element(tmp.begin(),tmp.end())<eps;