Select Git revision
correlation.h
correlation.h 17.05 KiB
#ifndef H_COMMON_WAVE_CORRELATION
#define H_COMMON_WAVE_CORRELATION
//Standard Library
#include <numeric>
//Gmsh Library
#include "gmsh.h"
//GmshFEM Library
#include "GmshFem.h"
#include "ExecutionTree.h"
#include "FieldInterface.h"
#include "FieldEvaluator.h"
#include "FieldNode.h"
#include "IndiceBucket.h"
#include "MathObject.h"
#include "OmpInterface.h"
#include "ScalarFunction.h"
#include "OperationsInterface.h"
#include "FieldForm0.h"
#include "Function.h"
#include "Domain.h"
#include "CSVio.h"
#include "MathObject.h"
//GmshFWI Library
#include "element.h"
/*
* (Optimized) Correlation function
*/
namespace gmshfem
{
namespace function
{
template<Physic T_Physic>
class BiWaveMultiFieldOperation : public OperationsInterface< std::complex<double>, Degree::Degree0 >
{
public:
typedef std::complex<double> D_Scalar;
constexpr static Degree D_Degree = Degree::Degree0;
BiWaveMultiFieldOperation() :
OperationsInterface< std::complex<double>, Degree::Degree0 >()
{
}
BiWaveMultiFieldOperation(const BiWaveMultiFieldOperation &other) :
OperationsInterface< std::complex<double>, Degree::Degree0 >(other)
{
}
virtual void operator()(const WaveMultiField<T_Physic> *fieldA, const WaveMultiField<T_Physic> *fieldB, std::vector< typename MathObject< std::complex<double>, Degree::Degree0 >::Object > &values, const std::vector< scalar::Precision< std::complex<double> > > &points, const std::vector< scalar::Precision< std::complex<double> > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const = 0;
};
template< class T_Op, Physic T_Physic >
class BiWaveMultiFieldNode final : public ExecutionTreeWithDegreeAndLeaves< typename T_Op::D_Scalar, T_Op::D_Degree >
{
protected:
T_Op _op;
const WaveMultiField<T_Physic> *_fieldA;
const WaveMultiField<T_Physic> *_fieldB;
public:
BiWaveMultiFieldNode(const WaveMultiField<T_Physic> *fieldA, const WaveMultiField<T_Physic> *fieldB) :
_op(), _fieldA(fieldA), _fieldB(fieldB)
{ }
BiWaveMultiFieldNode(const T_Op &other, const WaveMultiField<T_Physic> *fieldA, const WaveMultiField<T_Physic> *fieldB) :
_op(other), _fieldA(fieldA), _fieldB(fieldB)
{
}
virtual ~BiWaveMultiFieldNode() {}
virtual bool isConstant(const std::pair< int, int > &entity) const override
{
return false;
}
virtual void evaluate(std::vector< typename MathObject< typename T_Op::D_Scalar, T_Op::D_Degree >::Object > &values, const std::vector< scalar::Precision< typename T_Op::D_Scalar > > &points, const std::vector< scalar::Precision< typename T_Op::D_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const override
{
#pragma omp single
{
if(this->_memorize) {
this->_isInMemory = _op.getFromMemory(values, gaussPoints, elementType, entity);
}
}
if(!this->_isInMemory) {
#pragma omp single
values.resize(points.size() / 3);
_op(_fieldA, _fieldB, values, points, gaussPoints, elementType, entity);
#pragma omp barrier
if(this->_memorize) {
#pragma omp single
_op.setInMemory(values, gaussPoints, elementType, entity);
}
}
}
virtual BiWaveMultiFieldNode< T_Op,T_Physic > *copy() const override
{
BiWaveMultiFieldNode< T_Op,T_Physic > *cp = new BiWaveMultiFieldNode< T_Op,T_Physic >(_op, _fieldA, _fieldB);
if(this->_memorize) cp->activateMemorization();
return cp;
}
};
/* Template specialization should not appear here .., à voir avec Antho pour la généralisation ... */
template< Physic T_Physic >
class Correlate final : public BiWaveMultiFieldOperation< T_Physic >
{};
/* acoustic */
template<>
class Correlate<Physic::acoustic> final : public BiWaveMultiFieldOperation<Physic::acoustic>
{
private:
std::vector<unsigned int> _idx;
mutable std::vector< std::pair< int, int > > _keys;
mutable std::vector< std::complex<double> > _dofsValuesA;
mutable std::vector< std::complex<double> > _dofsValuesB;
mutable std::vector< scalar::Precision< std::complex<double> > > _functions;
mutable std::vector< int > _fsIndex;
mutable unsigned int _nbrDofsByElements;
mutable std::vector< double > _gmshJacobians, _gmshDeterminants, _gmshPoints, _gmshGaussPoints;
mutable std::vector< scalar::Precision< std::complex<double> > > _jacobians, _determinants;
mutable term::evaluator::FieldEvaluator< std::complex<double>, field::Form::Form0 > _fieldEvaluator;
public:
Correlate<Physic::acoustic>(const std::vector<unsigned int>& idx): _idx(idx) {}
Correlate<Physic::acoustic>(const Correlate &other) :
BiWaveMultiFieldOperation< Physic::acoustic >(other), _idx(other._idx) {
}
void operator()(const WaveMultiField<Physic::acoustic> *fieldA, const WaveMultiField<Physic::acoustic> *fieldB, std::vector< typename MathObject< std::complex<double>, Degree::Degree0 >::Object > &values, const std::vector< scalar::Precision< std::complex<double> > > &points, const std::vector< scalar::Precision< std::complex<double> > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const override
{
{
const unsigned int nbrOfElements = points.size() / gaussPoints.size();
const unsigned int nbrThreads = omp::nbrThreads();
const unsigned int myThreadID = omp::threadID();
#pragma omp single
_nbrDofsByElements = (*fieldA)[0].getFunctionSpace()->getBasisFunctions(false, _functions, _fsIndex, gaussPoints, elementType, entity);
std::vector< scalar::Precision< std::complex<double> > > coordinates;
//omp::setNested(true);
#pragma omp single
((*fieldA)[0].getFunctionSpace())->getKeys(false, _keys, coordinates, elementType, entity);
//omp::setNested(false);
const unsigned int begin = (myThreadID * _keys.size()) / nbrThreads;
const unsigned int end = ((myThreadID + 1) * _keys.size()) / nbrThreads;
const bool needOffset = !(*fieldA)[0].getFunctionSpace()->isConstantByElements();
#pragma omp single
_dofsValuesA.resize(_keys.size());
#pragma omp single
_dofsValuesB.resize(_keys.size());
const unsigned int nbrGaussPoints = gaussPoints.size() / 3;
#pragma omp single
scalar::copy(_gmshGaussPoints, gaussPoints);
if(_fieldEvaluator.needJacobians())
{
#pragma omp single
gmsh::model::mesh::preallocateJacobians(elementType, nbrGaussPoints, true, true, false, _gmshJacobians, _gmshDeterminants, _gmshPoints, entity.second);
gmsh::model::mesh::getJacobians(elementType, _gmshGaussPoints, _gmshJacobians, _gmshDeterminants, _gmshPoints, entity.second, myThreadID, nbrThreads);
#pragma omp barrier
#pragma omp single
scalar::move(_jacobians, _gmshJacobians);
#pragma omp single
scalar::move(_determinants, _gmshDeterminants);
}
for(auto it = _idx.begin(); it != _idx.end(); it++)
{
(*fieldA)[*it].getValues(_keys, _dofsValuesA, begin, end);
(*fieldB)[*it].getValues(_keys, _dofsValuesB, begin, end);
#pragma omp barrier
#pragma omp for
for(unsigned int i = 0; i < nbrOfElements; ++i)
{
Eigen::VectorX< std::complex<double> > vecDofA = Eigen::Map< const Eigen::VectorX< std::complex<double> > >(&(_dofsValuesA[i * _nbrDofsByElements]), _nbrDofsByElements);
Eigen::VectorX< std::complex<double> > vecDofB = Eigen::Map< const Eigen::VectorX< std::complex<double> > >(&(_dofsValuesB[i * _nbrDofsByElements]), _nbrDofsByElements);
for(unsigned int j = 0; j < nbrGaussPoints; ++j)
{
Eigen::MatrixX< scalar::Precision< std::complex<double>>> valueBasis = _fieldEvaluator(&_functions[((needOffset ? _fsIndex[i] * nbrGaussPoints : 0) + j) * field::NumberOfComponentsOfForm< field::Form::Form0 >::value * _nbrDofsByElements], &_jacobians[i * 9 * nbrGaussPoints + j * 9], &_determinants[i * nbrGaussPoints + j], _nbrDofsByElements);
values[i * nbrGaussPoints + j] += (valueBasis * vecDofA)(0) * (valueBasis * vecDofB)(0);
}
}
}
}
}
};
} // namespace function
} // namespace gmshfem
template<Physic T_Physic>ModelFunction correlate(const WaveMultiField<T_Physic>& u, const WaveMultiField<T_Physic>& v, const std::vector<unsigned int>& index)
{
return gmshfem::function::Function< std::complex<double>, gmshfem::Degree::Degree0 >( new gmshfem::function::BiWaveMultiFieldNode< gmshfem::function::Correlate< T_Physic >,T_Physic >( gmshfem::function::Correlate< T_Physic >(index), &u, &v) );
};
template<Physic T_Physic>
ModelFunction correlate(const WaveMultiField<T_Physic>& u, const WaveMultiField<T_Physic>& v)
{
std::vector<unsigned int> idx(u.size());
std::iota(idx.begin(), idx.end(), 0);
return correlate(u,v,idx);
};
/*
* (Optimized) Autocorrelate function
*/
namespace gmshfem
{
namespace function
{
template<Physic T_Physic>
class WaveMultiFieldOperation : public OperationsInterface< std::complex<double>, Degree::Degree0 >
{
public:
typedef std::complex<double> D_Scalar;
constexpr static Degree D_Degree = Degree::Degree0;
WaveMultiFieldOperation() :
OperationsInterface< std::complex<double>, Degree::Degree0 >()
{
}
WaveMultiFieldOperation(const WaveMultiFieldOperation &other) :
OperationsInterface< std::complex<double>, Degree::Degree0 >(other)
{
}
virtual void operator()(const WaveMultiField<T_Physic>* field, std::vector< typename MathObject< std::complex<double>, Degree::Degree0 >::Object > &values, const std::vector< scalar::Precision< std::complex<double> > > &points, const std::vector< scalar::Precision< std::complex<double> > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const = 0;
};
template< class T_Op, Physic T_Physic >
class WaveMultiFieldNode final : public ExecutionTreeWithDegreeAndLeaves< typename T_Op::D_Scalar, T_Op::D_Degree >
{
protected:
T_Op _op;
const WaveMultiField<T_Physic> *_field;
public:
WaveMultiFieldNode(const WaveMultiField<T_Physic>* field) :
_op(), _field(field)
{
}
WaveMultiFieldNode(const T_Op &other, const WaveMultiField<T_Physic>* field) :
_op(other), _field(field)
{
}
virtual ~WaveMultiFieldNode() {}
virtual bool isConstant(const std::pair< int, int > &entity) const override
{
return false;
}
virtual void evaluate(std::vector< typename MathObject< typename T_Op::D_Scalar, T_Op::D_Degree >::Object > &values, const std::vector< scalar::Precision< typename T_Op::D_Scalar > > &points, const std::vector< scalar::Precision< typename T_Op::D_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const override
{
#pragma omp single
{ if(this->_memorize) {
this->_isInMemory = _op.getFromMemory(values, gaussPoints, elementType, entity);
}
}
if(!this->_isInMemory) {
#pragma omp single
values.resize(points.size() / 3);
_op(_field, values, points, gaussPoints, elementType, entity);
#pragma omp barrier
if(this->_memorize) {
#pragma omp single
_op.setInMemory(values, gaussPoints, elementType, entity);
}
}
}
virtual WaveMultiFieldNode< T_Op,T_Physic > *copy() const override
{
WaveMultiFieldNode< T_Op,T_Physic > *cp = new WaveMultiFieldNode< T_Op,T_Physic >(_op, _field);
if(this->_memorize) cp->activateMemorization();
return cp;
}
};
/* Template specialization should not appear here .., à voir avec Antho pour la généralisation ... */
template< Physic T_Physic >
class Autocorrelate final : public WaveMultiFieldOperation< T_Physic >
{};
/* acoustic */
template<>
class Autocorrelate<Physic::acoustic> final : public WaveMultiFieldOperation<Physic::acoustic>
{
private:
std::vector<unsigned int> _idx;
mutable std::vector< std::pair< int, int > > _keys;
mutable std::vector< std::complex<double> > _dofsValues;
mutable std::vector< scalar::Precision< std::complex<double> > > _functions;
mutable std::vector< int > _fsIndex;
mutable unsigned int _nbrDofsByElements;
mutable std::vector< double > _gmshJacobians, _gmshDeterminants, _gmshPoints, _gmshGaussPoints;
mutable std::vector< scalar::Precision< std::complex<double> > > _jacobians, _determinants;
mutable term::evaluator::FieldEvaluator< std::complex<double>, field::Form::Form0 > _fieldEvaluator;
public:
Autocorrelate<Physic::acoustic>(const std::vector<unsigned int>& idx): _idx(idx)
{
}
Autocorrelate<Physic::acoustic>(const Autocorrelate &other) :
WaveMultiFieldOperation< Physic::acoustic >(other), _idx(other._idx)
{
}
void operator()(const WaveMultiField<Physic::acoustic>* field, std::vector< typename MathObject< std::complex<double>, Degree::Degree0 >::Object > &values, const std::vector< scalar::Precision< std::complex<double> > > &points, const std::vector< scalar::Precision< std::complex<double> > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const override
{
{
const unsigned int nbrOfElements = points.size() / gaussPoints.size();
const unsigned int nbrThreads = omp::nbrThreads();
const unsigned int myThreadID = omp::threadID();
#pragma omp single
_nbrDofsByElements = (*field)[0].getFunctionSpace()->getBasisFunctions(false, _functions, _fsIndex, gaussPoints, elementType, entity);
std::vector< scalar::Precision< std::complex<double> > > coordinates;
//omp::setNested(true);
#pragma omp single
(*field)[0].getFunctionSpace()->getKeys(false, _keys, coordinates, elementType, entity); //omp::setNested(false);
const unsigned int begin = (myThreadID * _keys.size()) / nbrThreads;
const unsigned int end = ((myThreadID + 1) * _keys.size()) / nbrThreads;
const bool needOffset = !(*field)[0].getFunctionSpace()->isConstantByElements();
#pragma omp single
_dofsValues.resize(_keys.size());
const unsigned int nbrGaussPoints = gaussPoints.size() / 3;
#pragma omp single
scalar::copy(_gmshGaussPoints, gaussPoints);
if(_fieldEvaluator.needJacobians())
{
#pragma omp single
gmsh::model::mesh::preallocateJacobians(elementType, nbrGaussPoints, true, true, false, _gmshJacobians, _gmshDeterminants, _gmshPoints, entity.second);
gmsh::model::mesh::getJacobians(elementType, _gmshGaussPoints, _gmshJacobians, _gmshDeterminants, _gmshPoints, entity.second, myThreadID, nbrThreads);
#pragma omp barrier
#pragma omp single
scalar::move(_jacobians, _gmshJacobians);
#pragma omp single
scalar::move(_determinants, _gmshDeterminants);
}
for(auto it = _idx.begin(); it != _idx.end(); it++)
{
(*field)[*it].getValues(_keys, _dofsValues, begin, end);
#pragma omp barrier
#pragma omp for
for(unsigned int i = 0; i < nbrOfElements; ++i)
{
Eigen::VectorX< std::complex<double> > vecDof = Eigen::Map< const Eigen::VectorX< std::complex<double> > >(&(_dofsValues[i * _nbrDofsByElements]), _nbrDofsByElements);
for(unsigned int j = 0; j < nbrGaussPoints; ++j)
{
Eigen::MatrixX< scalar::Precision< std::complex<double> > > valueBasis = _fieldEvaluator(&_functions[((needOffset ? _fsIndex[i] * nbrGaussPoints : 0) + j) * field::NumberOfComponentsOfForm< field::Form::Form0 >::value * _nbrDofsByElements], &_jacobians[i * 9 * nbrGaussPoints + j * 9], &_determinants[i * nbrGaussPoints + j], _nbrDofsByElements);
values[i * nbrGaussPoints + j] += (valueBasis * vecDof)(0) * std::conj((valueBasis * vecDof)(0));
}
}
}
}
}
};
} // namespace function
} // namespace gmshfem
template<Physic T_Physic>
ModelFunction autocorrelate(const WaveMultiField<T_Physic>& u, const std::vector<unsigned int>& index)
{
return gmshfem::function::Function< std::complex<double>, gmshfem::Degree::Degree0 >( new gmshfem::function::WaveMultiFieldNode< gmshfem::function::Autocorrelate< T_Physic >,T_Physic >( gmshfem::function::Autocorrelate< T_Physic >(index), &u) );
};
template<Physic T_Physic>
ModelFunction autocorrelate(const WaveMultiField<T_Physic>& u)
{
std::vector<unsigned int> idx(u.size());
std::iota(idx.begin(), idx.end(), 0);
return autocorrelate(u,idx);
};
#endif // H_COMMON_WAVE_CORRELATION