Skip to content
Snippets Groups Projects
Commit 9749c58e authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

dgAlgorithm compilation ok

parent fee0b8bc
No related branches found
No related tags found
No related merge requests found
......@@ -60,7 +60,7 @@ void fullMatrix<std::complex<double> >::mult(const fullMatrix<std::complex<doubl
}
template<>
void fullMatrix<double>::gemm(fullMatrix<double> &a, fullMatrix<double> &b,
void fullMatrix<double>::gemm(const fullMatrix<double> &a, const fullMatrix<double> &b,
double alpha, double beta)
{
int M = size1(), N = size2(), K = a.size2();
......@@ -70,8 +70,8 @@ void fullMatrix<double>::gemm(fullMatrix<double> &a, fullMatrix<double> &b,
}
template<>
void fullMatrix<std::complex<double> >::gemm(fullMatrix<std::complex<double> > &a,
fullMatrix<std::complex<double> > &b,
void fullMatrix<std::complex<double> >::gemm(const fullMatrix<std::complex<double> > &a,
const fullMatrix<std::complex<double> > &b,
std::complex<double> alpha,
std::complex<double> beta)
{
......
......@@ -143,7 +143,7 @@ class fullMatrix
}
#endif
;
void gemm(fullMatrix<scalar> &a, fullMatrix<scalar> &b,
void gemm(const fullMatrix<scalar> &a, const fullMatrix<scalar> &b,
scalar alpha=1., scalar beta=1.)
#if !defined(HAVE_BLAS)
{
......
......@@ -11,6 +11,7 @@ set(SRC
SElement.cpp
eigenSolver.cpp
dgGroupOfElements.cpp
dgAlgorithm.cpp
)
file(GLOB HDR RELATIVE ${CMAKE_SOURCE_DIR}/Solver *.h)
......
#include "dgAlgorithm.h"
#include "dgGroupOfElements.h"
#include "dgConservationLaw.h"
/*
compute
\int \vec{f} \cdot \grad \phi dv
*/
void dgAlgorithm::residualVolume ( dofManager &dof, // the DOF manager (maybe useless here)
void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
const groupOfElements & group ) // the residual
const dgGroupOfElements & group ) // the residual
{
// ----- 1 ---- get the solution at quadrature points
// ----- 1.1 --- allocate a matrix of size (nbFields * nbElements, nbQuadraturePoints)
fullMatrix<double> solutionQP (group.getNbElements() * group.getNbFields(), group.getNbIntegrationPoints());
int nbFields = group.getNbFields();
fullMatrix<double> solutionQP (group.getNbElements() * nbFields, group.getNbIntegrationPoints());
// ----- 1.2 --- multiply the solution by the collocation matrix
group.getSolution().mult ( group.getCollocationMatrix() , solutionQP);
// ----- 1.3 --- if the conservation law is diffusive, compute the gradients too
fullMatrix<double> *gradientSolutionQP= 0;
if (claw.diffusiveFlux()){
gradientSolutionQP = new fullMatrix<double> (group.getNbElements() * group.getNbFields() * 3, group.getNbIntegrationPoints());
gradientSolutionQP = new fullMatrix<double> (group.getNbElements() * nbFields * 3, group.getNbIntegrationPoints());
group.getGradientOfSolution().mult ( group.getCollocationMatrix() , *gradientSolutionQP);
}
// ----- 2 ---- compute all fluxes (diffusive and convective) at integration points
// ----- 2.1 --- allocate elementary fluxes (compued in XYZ coordinates)
fullMatrix<double> fConv[3] = {fullMatrix<double>( group.getNbIntegrationPoints(), group.nbFields() ),
fullMatrix<double>( group.getNbIntegrationPoints(), group.nbFields() ),
fullMatrix<double>( group.getNbIntegrationPoints(), group.nbFields() )};
fullMatrix<double> fDiff[3] = {fullMatrix<double>( group.getNbIntegrationPoints(), group.nbFields() ),
fullMatrix<double>( group.getNbIntegrationPoints(), group.nbFields() ),
fullMatrix<double>( group.getNbIntegrationPoints(), group.nbFields() )};
fullMatrix<double> fConv[3] = {fullMatrix<double>( group.getNbIntegrationPoints(), nbFields ),
fullMatrix<double>( group.getNbIntegrationPoints(), nbFields ),
fullMatrix<double>( group.getNbIntegrationPoints(), nbFields )};
fullMatrix<double> fDiff[3] = {fullMatrix<double>( group.getNbIntegrationPoints(), nbFields ),
fullMatrix<double>( group.getNbIntegrationPoints(), nbFields ),
fullMatrix<double>( group.getNbIntegrationPoints(), nbFields )};
// ----- 2.2 --- allocate parametric fluxes (compued in UVW coordinates) for all elements at all integration points
fullMatrix<double> Fuvw[3] = {fullMatrix<double> (group.getNbElements() * group.getNbFields(), group.getNbIntegrationPoints()),
fullMatrix<double> (group.getNbElements() * group.getNbFields(), group.getNbIntegrationPoints()),
fullMatrix<double> (group.getNbElements() * group.getNbFields(), group.getNbIntegrationPoints())};
fullMatrix<double> Fuvw[3] = {fullMatrix<double> (group.getNbElements() * nbFields, group.getNbIntegrationPoints()),
fullMatrix<double> (group.getNbElements() * nbFields, group.getNbIntegrationPoints()),
fullMatrix<double> (group.getNbElements() * nbFields, group.getNbIntegrationPoints())};
// ----- 2.3 --- iterate on elements
for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
// ----- 2.3.1 --- build a small object that contains elementary solution, jacobians, gmsh element
fullMatrix<double> solutionQPe (solutionQP, iElement*claw.nbFields(),claw.nbFields() );
fullMatrix<double> *gradSolutionQPe;
if (claw.diffusiveFlux()) gradSolutionQPe = new fullMatrix<double>(*gradSolutionQP, 3*iElement*claw.nbFields(),3*claw.nbFields() );
if (claw.diffusiveFlux()) gradSolutionQPe = new fullMatrix<double>(*gradientSolutionQP, 3*iElement*claw.nbFields(),3*claw.nbFields() );
else gradSolutionQPe = new fullMatrix<double>;
dgElement DGE( group.getElement(iElement), solutionQPe, *gradSolutionQPe, group.getIntegrationPointsMatrix());
// ----- 2.3.2 --- compute fluxes in XYZ coordinates
......@@ -59,7 +62,7 @@ void dgAlgorithm::residualVolume ( dofManager &dof, // the DOF manager (maybe us
const double detJ = group.getDetJ (iElement, iPt);
const double factor = invJ * detJ;
// compute fluxes in the reference coordinate system at this integration point
for (int k=0;k<group.nbFields();k++) {
for (int k=0;k<nbFields;k++) {
fuvwe(iPt,k) += ( fConv[iXYZ](iPt,k) + fDiff[iXYZ](iPt,k)) * factor;
}
}
......@@ -68,5 +71,5 @@ void dgAlgorithm::residualVolume ( dofManager &dof, // the DOF manager (maybe us
}
// ----- 3 ---- do the redistribution at nodes using as many BLAS3 operations as there are local coordinates
for (int iUVW=0;iUVW<group.getDimUVW();iUVW++)
group.getResidual().dgemm(group.getRedistributionMatrix(iUVW), Fuvw[iUVW]);
group.getResidual().gemm(group.getFluxRedistributionMatrix(iUVW), Fuvw[iUVW]);
}
......@@ -3,27 +3,21 @@
#include "dofManager.h"
class groupOfElements;
class groupOfFaces;
class dgGroupOfElements;
class dgGroupOfFaces;
class dgConservationLaw;
class dgAlgorithm () {
class dgAlgorithm {
public :
void residualVolume ( dofManager &dof,
dgConservationLaw &law,
groupOfElements & group ,
fullMatrix<double> &solution,
fullMatrix<double> &residual);
void residualInterface ( dofManager &dof,
dgConservationLaw &law,
groupOfFaces & group ,
fullMatrix<double> &solution,
fullMatrix<double> &residual);
void residualBoundaryConditions ( dofManager &dof,
dgConservationLaw &law,
groupOfElements & group ,
fullMatrix<double> &solution,
fullMatrix<double> &residual);
void residualVolume ( /*dofManager &dof,*/
const dgConservationLaw &law,
const dgGroupOfElements & group);
void residualInterface ( /*dofManager &dof,*/
const dgConservationLaw &law,
const dgGroupOfFaces & group);
void residualBoundaryConditions ( /*dofManager &dof,*/
const dgConservationLaw &law,
const dgGroupOfElements & group);
};
......
......@@ -31,13 +31,13 @@ public:
_riemannSolver(0),_maxConvectiveSpeed (0) {}
~dgConservationLaw () {}
int nbFields() const {return nbf;}
int nbFields() const {return _nbf;}
inline const dgTerm * convectiveFlux () {return _convective;}
inline const dgTerm * diffusiveFlux () {return _diffusive;}
inline const dgTerm * sourceTerm () {return _source;}
inline const dgFaceTerm * riemannSolver () {return _riemannSolver;}
inline const dgTerm * maxConvectiveSpeed () {return _maxConvectiveSpeed;}
inline const dgTerm * convectiveFlux () const {return _convective;}
inline const dgTerm * diffusiveFlux () const {return _diffusive;}
inline const dgTerm * sourceTerm () const {return _source;}
inline const dgFaceTerm * riemannSolver () const {return _riemannSolver;}
inline const dgTerm * maxConvectiveSpeed () const {return _maxConvectiveSpeed;}
};
......
......@@ -19,14 +19,14 @@ class MEdge;
class functionSpace;
class dgElement {
MElement *_element;
const MElement *_element;
// solution at points
const fullMatrix<double> &_solution, &_integration, &_gradients;
public:
dgElement (MElement *e, const fullMatrix<double> &sol, const fullMatrix<double> &integ)
dgElement (const MElement *e, const fullMatrix<double> &sol, const fullMatrix<double> &integ)
: _element(e), _solution(sol), _integration(integ), _gradients(sol)
{}
dgElement (MElement *e, const fullMatrix<double> &sol, const fullMatrix<double> &grads, const fullMatrix<double> &integ)
dgElement (const MElement *e, const fullMatrix<double> &sol, const fullMatrix<double> &grads, const fullMatrix<double> &integ)
: _element(e), _solution(sol), _integration(integ), _gradients(grads)
{}
};
......@@ -76,11 +76,11 @@ public:
inline const fullMatrix<double> & getCollocationMatrix () const {return *_collocation;}
inline const fullMatrix<double> & getFluxRedistributionMatrix (int i) const {return *_redistributionFluxes[i];}
inline const fullMatrix<double> & getSourceRedistributionMatrix () const {return *_redistributionSource;}
inline const fullMatrix<double> & getSolution () const {return *_solution;}
inline const fullMatrix<double> & getGradientOfSolution () const {return *_gradSolution;}
inline fullMatrix<double> & getSolution () const {return *_solution;}
inline fullMatrix<double> & getGradientOfSolution () const {return *_gradSolution;}
// get a proxy on the solution for element iElement
inline fullMatrix<double> getSolution (int iElement) const {return fullMatrix<double>(*_solution, iElement*_nbFields, _nbFields);}
inline const fullMatrix<double> & getResidual () const {return *_solution;}
inline fullMatrix<double> & getResidual () const {return *_solution;}
// get a proxy on the residual for element iElement
inline fullMatrix<double> getResidual (int iElement) const {return fullMatrix<double>(*_residual, iElement*_nbFields, _nbFields);}
inline double getDetJ (int iElement, int iGaussPoint) const {return (*_mapping)(iElement, 10*iGaussPoint + 9);}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment