1 // Copyright (C) 2007-2023 CEA, EDF
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 (EDF R&D)
21 #include "MEDCouplingFieldDiscretizationOnNodesFE.hxx"
22 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
23 #include "InterpKernelDenseMatrix.hxx"
24 #include "InterpKernelRootsMultiDim.hxx"
25 #include "MEDCouplingUMesh.hxx"
26 #include "InterpolationHelper.txx"
27 #include "InterpKernelGaussCoords.hxx"
31 using namespace MEDCoupling;
33 const char MEDCouplingFieldDiscretizationOnNodesFE::REPR[]="FE";
35 std::string MEDCouplingFieldDiscretizationOnNodesFE::getStringRepr() const
37 return std::string(REPR);
40 void MEDCouplingFieldDiscretizationOnNodesFE::reprQuickOverview(std::ostream& stream) const
42 stream << "NodeFE spatial discretization.";
45 MCAuto<MEDCouplingFieldDiscretization> MEDCouplingFieldDiscretizationOnNodesFE::aggregate(std::vector<const MEDCouplingFieldDiscretization *>& fds) const
47 return EasyAggregate<MEDCouplingFieldDiscretizationOnNodesFE>(fds);
50 bool MEDCouplingFieldDiscretizationOnNodesFE::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
54 reason="other spatial discretization is NULL, and this spatial discretization (Node FE) is defined.";
57 const MEDCouplingFieldDiscretizationOnNodesFE *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationOnNodesFE *>(other);
60 reason="Spatial discrtization of this is ON_NODES_FE, which is not the case of other.";
65 * This method is simply called by MEDCouplingFieldDiscretization::deepCopy. It performs the deep copy of \a this.
67 * \sa MEDCouplingFieldDiscretization::deepCopy.
69 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationOnNodesFE::clone() const
71 return new MEDCouplingFieldDiscretizationOnNodesFE;
74 void MEDCouplingFieldDiscretizationOnNodesFE::checkCompatibilityWithNature(NatureOfField nat) const
76 if(nat!=IntensiveMaximum)
77 throw INTERP_KERNEL::Exception("Invalid nature for NodeFE field : expected IntensiveMaximum !");
80 MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationOnNodesFE::getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const
83 throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationOnNodesFE::getMeasureField : mesh instance specified is NULL !");
84 throw INTERP_KERNEL::Exception("getMeasureField on MEDCouplingFieldDiscretizationOnNodesFE : not implemented yet !");
90 const MEDCouplingGaussLocalization* _gl;
91 std::size_t _nb_pts_in_cell;
92 const double *_pts_in_cell;
95 Functor(const MEDCouplingGaussLocalization& gl, std::size_t nbPtsInCell, const double *ptsInCell, const double point[3]):_gl(&gl),_nb_pts_in_cell(nbPtsInCell),
96 _pts_in_cell(ptsInCell),_point(point) { }
97 std::vector<double> operator()(const std::vector<double>& x)
99 MEDCouplingGaussLocalization gl(_gl->getType(),_gl->getRefCoords(),x,{1.0});
100 MCAuto<DataArrayDouble> shapeFunc = gl.getShapeFunctionValues();
101 const double *shapeFuncPtr( shapeFunc->begin() );
102 std::vector<double> ret(3,0);
103 for(std::size_t iPt = 0; iPt < _nb_pts_in_cell ; ++iPt)
105 for(short iDim = 0 ; iDim < 3 ; ++iDim)
106 ret[iDim] += shapeFuncPtr[iPt] * _pts_in_cell[3*iPt + iDim];
108 ret[0] -= _point[0]; ret[1] -= _point[1]; ret[2] -= _point[2];
113 bool IsInside3D(const MEDCouplingGaussLocalization& gl, const std::vector<double>& ptsInCell, const double locInReal[3], double locInRef[3])
115 constexpr double EPS_IN_OUT = 1e-12;
116 std::size_t nbPtsInCell(ptsInCell.size()/3);
118 const double *refCoo(gl.getRefCoords().data());
119 INTERP_KERNEL::NormalizedCellType ct(gl.getType());
120 Functor func(gl,nbPtsInCell,ptsInCell.data(),locInReal);
122 auto myJacobian = [&gl,nbPtsInCell,ptsInCell](const std::vector<double>& x, const std::vector<double>&, INTERP_KERNEL::DenseMatrix& jacobian)
124 MEDCouplingGaussLocalization mygl(gl.getType(),gl.getRefCoords(),x,{1.0});
125 MCAuto<DataArrayDouble> shapeFunc = mygl.getDerivativeOfShapeFunctionValues();
126 for(std::size_t i = 0 ; i < 3 ; ++i)
127 for(std::size_t j = 0 ; j < 3 ; ++j)
130 for( std::size_t k = 0 ; k < nbPtsInCell ; ++k )
131 res += ptsInCell[k*3+i] * shapeFunc->getIJ(0,3*k+j);
132 jacobian[ i ][ j ] = res;
136 // loop on refcoords as initialization point for Newton algo. vini is the initialization vector of Newton.
137 for(std::size_t attemptId = 0 ; attemptId < nbPtsInCell ; ++attemptId)
139 std::vector<double> vini(refCoo + attemptId*3, refCoo + (attemptId+1)*3);
143 //INTERP_KERNEL::SolveWithNewton(vini,check,func);
144 INTERP_KERNEL::SolveWithNewtonWithJacobian(vini,check,func,myJacobian);
145 ret = (check==false);//looks strange but OK regarding newt (SolveWithNewton) at page 387 of numerical recipes for semantic of check parameter
147 catch( INTERP_KERNEL::Exception& ex )
148 { ret = false; }// Something get wrong during Newton process
150 {//Newton has converged. Now check if it converged to a point inside cell
151 if( ! INTERP_KERNEL::GaussInfo::IsInOrOutForReference(ct,vini.data(),EPS_IN_OUT) )
152 {// converged but locInReal has been detected outside of cell
158 locInRef[0] = vini[0]; locInRef[1] = vini[1]; locInRef[2] = vini[2];
162 std::fill(locInRef,locInRef+3,std::numeric_limits<double>::max());
166 void MEDCouplingFieldDiscretizationOnNodesFE::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *ptsCoo, double *res) const
168 MCAuto<DataArrayDouble> res2( this->getValueOnMulti(arr,mesh,ptsCoo,1) );
169 std::copy(res2->begin(),res2->end(),res);
172 void MEDCouplingFieldDiscretizationOnNodesFE::GetRefCoordOfListOf3DPtsIn3D(const MEDCouplingUMesh *umesh, const double *ptsCoo, mcIdType nbOfPts,
173 std::function<void(const MEDCouplingGaussLocalization&, const std::vector<mcIdType>&)> customFunc)
175 const double *coordsOfMesh( umesh->getCoords()->begin() );
176 MEDCouplingNormalizedUnstructuredMesh<3,3> mesh_wrapper(umesh);
177 BBTreeStandAlone<3,mcIdType> tree( INTERP_KERNEL::BuildBBTree( mesh_wrapper ) );
178 for(mcIdType iPt = 0 ; iPt < nbOfPts ; ++iPt)
180 std::vector<mcIdType> elems;
181 tree.getElementsAroundPoint(ptsCoo+3*iPt,elems);
183 for(auto cellId = elems.cbegin() ; cellId != elems.cend() && !found ; ++cellId)
185 INTERP_KERNEL::NormalizedCellType gt( umesh->getTypeOfCell(*cellId) );
186 std::vector<mcIdType> conn;
187 umesh->getNodeIdsOfCell(*cellId,conn);
188 MCAuto<DataArrayDouble> refCoo( MEDCouplingGaussLocalization::GetDefaultReferenceCoordinatesOf(gt) );
189 std::vector<double> refCooCpp(refCoo->begin(),refCoo->end());
190 std::vector<double> gsCoo(ptsCoo + iPt*3,ptsCoo + (iPt+1)*3);
191 MEDCouplingGaussLocalization gl(gt,refCooCpp,{0,0,0},{1.});
192 std::vector<double> ptsInCell; ptsInCell.reserve(conn.size()*gl.getDimension());
193 std::for_each( conn.cbegin(), conn.cend(), [coordsOfMesh,&ptsInCell](mcIdType c) { ptsInCell.insert(ptsInCell.end(),coordsOfMesh+c*3,coordsOfMesh+(c+1)*3); } );
194 std::vector<double> locInRef(3);
195 if( IsInside3D(gl,ptsInCell,gsCoo.data(),locInRef.data()) )
197 gl.setGaussCoords(locInRef);
203 THROW_IK_EXCEPTION("getValueOnMulti on MEDCouplingFieldDiscretizationOnNodesFE : fail to locate point #" << iPt << " X=" << ptsCoo[0] << " Y=" << ptsCoo[1] << " Z=" << ptsCoo[2] << " !");
207 const MEDCouplingUMesh *MEDCouplingFieldDiscretizationOnNodesFE::checkConfig3D(const MEDCouplingMesh *mesh) const
209 const MEDCouplingUMesh *umesh( dynamic_cast<const MEDCouplingUMesh *>(mesh) );
211 THROW_IK_EXCEPTION("getValueOn : not implemented yet for type != MEDCouplingUMesh !");
212 if(umesh->getSpaceDimension() != 3 || umesh->getMeshDimension() != 3)
213 THROW_IK_EXCEPTION("getValueOn : implemented only for meshes with spacedim == 3 and meshdim == 3 !");
214 umesh->checkConsistency();
218 DataArrayDouble *MEDCouplingFieldDiscretizationOnNodesFE::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfTargetPoints) const
220 if(!arr || !arr->isAllocated())
221 throw INTERP_KERNEL::Exception("getValueOnMulti : input array is null or not allocated !");
222 mcIdType nbOfRows=getNumberOfMeshPlaces(mesh);
223 if(arr->getNumberOfTuples()!=nbOfRows)
225 THROW_IK_EXCEPTION( "getValueOnMulti : input array does not have correct number of tuples ! Excepted " << nbOfRows << " having " << arr->getNumberOfTuples() << " !")
227 const MEDCouplingUMesh *umesh = checkConfig3D(mesh);
228 std::size_t nbCompo( arr->getNumberOfComponents() );
229 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
230 ret->alloc(nbOfTargetPoints,nbCompo);
231 double *res( ret->getPointer() );
233 auto arrayFeeder = [arr,&res,nbCompo](const MEDCouplingGaussLocalization& gl, const std::vector<mcIdType>& conn)
235 MCAuto<DataArrayDouble> resVector( gl.getShapeFunctionValues() );
237 std::for_each(res,res+nbCompo,[](double& v) { v = 0.0; });
238 for(std::size_t iComp = 0 ; iComp < nbCompo ; ++iComp)
239 for(int iPt = 0 ; iPt < gl.getNumberOfPtsInRefCell(); ++iPt)
242 res[iComp] += resVector->getIJ(0,iPt) * arr->getIJ(conn[iPt],iComp);
249 GetRefCoordOfListOf3DPtsIn3D(umesh,loc,nbOfTargetPoints,arrayFeeder);
254 * Returns for each \a nbOfPoints point in \a loc its coordinate in reference element.
256 MCAuto<DataArrayDouble> MEDCouplingFieldDiscretizationOnNodesFE::getCooInRefElement(const MEDCouplingMesh *mesh, const double *loc, mcIdType nbOfPoints) const
258 const MEDCouplingUMesh *umesh = checkConfig3D(mesh);
259 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
260 ret->alloc(nbOfPoints,3);
261 double *retPtr(ret->getPointer() );
263 auto arrayFeeder = [&retPtr](const MEDCouplingGaussLocalization& gl, const std::vector<mcIdType>& conn)
265 std::vector<double> resVector( gl.getGaussCoords() );
267 std::copy(resVector.begin(),resVector.end(),retPtr);
272 GetRefCoordOfListOf3DPtsIn3D(umesh,loc,nbOfPoints,arrayFeeder);