-// Copyright (C) 2007-2022 CEA/DEN, EDF R&D
+// Copyright (C) 2007-2024 CEA, EDF
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// Author : Anthony Geay (EDF R&D)
#include "MEDCouplingFieldDiscretization.hxx"
+#include "MEDCouplingFieldDiscretizationOnNodesFE.hxx"
#include "MEDCouplingCMesh.hxx"
#include "MEDCouplingUMesh.hxx"
#include "MEDCouplingFieldDouble.hxx"
#include "InterpKernelAutoPtr.hxx"
#include "InterpKernelGaussCoords.hxx"
#include "InterpKernelMatrixTools.hxx"
+#include "InterpKernelDenseMatrix.hxx"
#include <set>
#include <list>
const char MEDCouplingFieldDiscretizationP0::REPR[]="P0";
-const TypeOfField MEDCouplingFieldDiscretizationP0::TYPE=ON_CELLS;
-
const char MEDCouplingFieldDiscretizationP1::REPR[]="P1";
-const TypeOfField MEDCouplingFieldDiscretizationP1::TYPE=ON_NODES;
-
const mcIdType MEDCouplingFieldDiscretizationPerCell::DFT_INVALID_LOCID_VALUE=-1;
const char MEDCouplingFieldDiscretizationGauss::REPR[]="GAUSS";
-const TypeOfField MEDCouplingFieldDiscretizationGauss::TYPE=ON_GAUSS_PT;
-
const char MEDCouplingFieldDiscretizationGaussNE::REPR[]="GSSNE";
-const TypeOfField MEDCouplingFieldDiscretizationGaussNE::TYPE=ON_GAUSS_NE;
-
const char MEDCouplingFieldDiscretizationKriging::REPR[]="KRIGING";
-const TypeOfField MEDCouplingFieldDiscretizationKriging::TYPE=ON_NODES_KR;
-
// doc is here http://www.code-aster.org/V2/doc/default/fr/man_r/r3/r3.01.01.pdf
const double MEDCouplingFieldDiscretizationGaussNE::FGP_POINT1[1]={0.};
const double MEDCouplingFieldDiscretizationGaussNE::FGP_SEG2[2]={1.,1.};
return new MEDCouplingFieldDiscretizationGaussNE;
case MEDCouplingFieldDiscretizationKriging::TYPE:
return new MEDCouplingFieldDiscretizationKriging;
+ case MEDCouplingFieldDiscretizationOnNodesFE::TYPE:
+ return new MEDCouplingFieldDiscretizationOnNodesFE;
default:
throw INTERP_KERNEL::Exception("Chosen discretization is not implemented yet.");
}
return MEDCouplingFieldDiscretizationGaussNE::TYPE;
if(repr==MEDCouplingFieldDiscretizationKriging::REPR)
return MEDCouplingFieldDiscretizationKriging::TYPE;
+ if(repr==MEDCouplingFieldDiscretizationOnNodesFE::REPR)
+ return MEDCouplingFieldDiscretizationOnNodesFE::TYPE;
throw INTERP_KERNEL::Exception("Representation does not match with any field discretization !");
}
std::string MEDCouplingFieldDiscretization::GetTypeOfFieldRepr(TypeOfField type)
{
- if(type==MEDCouplingFieldDiscretizationP0::TYPE)
- return MEDCouplingFieldDiscretizationP0::REPR;
- if(type==MEDCouplingFieldDiscretizationP1::TYPE)
- return MEDCouplingFieldDiscretizationP1::REPR;
- if(type==MEDCouplingFieldDiscretizationGauss::TYPE)
- return MEDCouplingFieldDiscretizationGauss::REPR;
- if(type==MEDCouplingFieldDiscretizationGaussNE::TYPE)
- return MEDCouplingFieldDiscretizationGaussNE::REPR;
- if(type==MEDCouplingFieldDiscretizationKriging::TYPE)
- return MEDCouplingFieldDiscretizationKriging::REPR;
- throw INTERP_KERNEL::Exception("GetTypeOfFieldRepr : Representation does not match with any field discretization !");
+ switch(type)
+ {
+ case MEDCouplingFieldDiscretizationP0::TYPE:
+ return MEDCouplingFieldDiscretizationP0::REPR;
+ case MEDCouplingFieldDiscretizationP1::TYPE:
+ return MEDCouplingFieldDiscretizationP1::REPR;
+ case MEDCouplingFieldDiscretizationGauss::TYPE:
+ return MEDCouplingFieldDiscretizationGauss::REPR;
+ case MEDCouplingFieldDiscretizationGaussNE::TYPE:
+ return MEDCouplingFieldDiscretizationGaussNE::REPR;
+ case MEDCouplingFieldDiscretizationKriging::TYPE:
+ return MEDCouplingFieldDiscretizationKriging::REPR;
+ case MEDCouplingFieldDiscretizationOnNodesFE::TYPE:
+ return MEDCouplingFieldDiscretizationOnNodesFE::REPR;
+ default:
+ throw INTERP_KERNEL::Exception("GetTypeOfFieldRepr : Representation does not match with any field discretization !");
+ }
}
bool MEDCouplingFieldDiscretization::isEqual(const MEDCouplingFieldDiscretization *other, double eps) const
}
}
-template<class FIELD_DISC>
-MCAuto<MEDCouplingFieldDiscretization> MEDCouplingFieldDiscretization::EasyAggregate(std::vector<const MEDCouplingFieldDiscretization *>& fds)
-{
- if(fds.empty())
- throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretization::aggregate : input array is empty");
- for(const MEDCouplingFieldDiscretization * it : fds)
- {
- const FIELD_DISC *itc(dynamic_cast<const FIELD_DISC *>(it));
- if(!itc)
- throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretization::aggregate : same field discretization expected for all input discretizations !");
- }
- return fds[0]->clone();
-}
-
MEDCouplingFieldDiscretization::~MEDCouplingFieldDiscretization()
{
}
{
if(!mesh)
throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getMeasureField : mesh instance specified is NULL !");
- MCAuto<MEDCouplingFieldDouble> vol=mesh->getMeasureField(isAbs);
- const double *volPtr=vol->getArray()->begin();
+ MCAuto<MEDCouplingUMesh> umesh(mesh->buildUnstructured());
+ const double *coordsOfMesh( umesh->getCoords()->begin() );
+ auto spaceDim(mesh->getSpaceDimension());
+ auto meshDim(mesh->getMeshDimension());
MCAuto<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(ON_GAUSS_PT);
ret->setMesh(mesh);
ret->setDiscretization(const_cast<MEDCouplingFieldDiscretizationGauss *>(this));
_discr_per_cell->checkAllocated();
if(_discr_per_cell->getNumberOfComponents()!=1)
throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getMeasureField : no discr per cell array defined but with nb of components different from 1 !");
- if(_discr_per_cell->getNumberOfTuples()!=vol->getNumberOfTuples())
+ if(_discr_per_cell->getNumberOfTuples()!=mesh->getNumberOfCells())
throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getMeasureField : no discr per cell array defined but mismatch between nb of cells of mesh and size of spatial disr array !");
MCAuto<DataArrayIdType> offset=getOffsetArr(mesh);
MCAuto<DataArrayDouble> arr=DataArrayDouble::New(); arr->alloc(getNumberOfTuples(mesh),1);
const MEDCouplingGaussLocalization& loc=_loc[locId];
mcIdType nbOfGaussPt=loc.getNumberOfGaussPt();
INTERP_KERNEL::AutoPtr<double> weights=new double[nbOfGaussPt];
- double sum=std::accumulate(loc.getWeights().begin(),loc.getWeights().end(),0.);
- std::transform(loc.getWeights().begin(),loc.getWeights().end(),(double *)weights,std::bind(std::multiplies<double>(),std::placeholders::_1,1./sum));
for(const mcIdType *cellId=curIds->begin();cellId!=curIds->end();cellId++)
- for(mcIdType j=0;j<nbOfGaussPt;j++)
- arrPtr[offsetPtr[*cellId]+j]=weights[j]*volPtr[*cellId];
+ {
+ std::vector<mcIdType> conn;
+ umesh->getNodeIdsOfCell(*cellId,conn);
+ std::vector<double> ptsInCell; ptsInCell.reserve(conn.size()*loc.getDimension());
+ std::for_each( conn.cbegin(), conn.cend(), [spaceDim,coordsOfMesh,&ptsInCell](mcIdType c) { ptsInCell.insert(ptsInCell.end(),coordsOfMesh+c*spaceDim,coordsOfMesh+(c+1)*spaceDim); } );
+ std::size_t nbPtsInCell(ptsInCell.size()/spaceDim);
+ INTERP_KERNEL::DenseMatrix jacobian(spaceDim,meshDim);
+ MCAuto<DataArrayDouble> shapeFunc = loc.getDerivativeOfShapeFunctionValues();
+ for(mcIdType iGPt = 0 ; iGPt < nbOfGaussPt ; ++iGPt)
+ {
+ for(auto i = 0 ; i < spaceDim ; ++i)
+ for(auto j = 0 ; j < meshDim ; ++j)
+ {
+ double res = 0.0;
+ for( std::size_t k = 0 ; k < nbPtsInCell ; ++k )
+ res += ptsInCell[spaceDim*k+i] * shapeFunc->getIJ(iGPt,meshDim*k+j);
+ jacobian[ i ][ j ] = res;
+ }
+ arrPtr[offsetPtr[*cellId]+iGPt]=std::abs( jacobian.toJacobian() )*loc.getWeight(FromIdType<int>(iGPt));
+ }
+ }
}
else
{
