[Transport] Write tensors in form of a dim*dim vector

* On release 2.6, dune-grid will throw an exception due to a wrong setup of the VTKInfo. So this will fail right now.
* A way to solve this would be to create a simple writer and add data together with the correct VTKInfo.
parent fe66a75f
......@@ -7,114 +7,24 @@
namespace Dune{
namespace Dorie{
// template<class T, class Parameter, class GFWaterFlux, class GFWaterContent, bool tensor_glyph = true>
// class EffectiveHydrodynamicDispersionAdapter
// : public Dune::PDELab::GridFunctionBase<
// Dune::PDELab::GridFunctionTraits<
// typename T::GridView,
// typename T::RangeField,
// tensor_glyph ? 3 : T::GridView::dimension,
// std::conditional_t<tensor_glyph,Dune::FieldMatrix<typename T::RangeField,3,3>,typename T::Tensor>
// >
// ,EffectiveHydrodynamicDispersionAdapter<T,Parameter,GFWaterFlux,GFWaterContent,tensor_glyph>
// >
// {
// public:
// using Traits = Dune::PDELab::GridFunctionTraits<
// typename T::GridView,
// typename T::RangeField,
// tensor_glyph ? 3 : T::GridView::dimension,
// std::conditional_t<tensor_glyph,Dune::FieldMatrix<typename T::RangeField,3,3>,typename T::Tensor>
// >;
// private:
// using RF = typename T::RangeField;
// using Domain = typename T::Domain;
// public:
// /*-------------------------------------------------------------------*//**
// * @brief Constructor for the EffectiveHydrodynamicDispersionAdapter
// *
// * @param gv GridView
// * @param p Parametrization class
// * @see RichardsEquationParameter
// */
// EffectiveHydrodynamicDispersionAdapter(const typename Traits::GridViewType& gv,
// std::shared_ptr<const Parameter> param,
// std::shared_ptr<const GFWaterFlux> gf_water_flux,
// std::shared_ptr<const GFWaterContent> gf_water_content)
// : _gv(gv)
// , _param(param)
// , _gf_water_flux(gf_water_flux)
// , _gf_water_content(gf_water_content)
// {}
// /*-------------------------------------------------------------------*//**
// * @brief Evaluation of the conductivity for a given entity e in an
// * entity position x
// *
// * @param[in] e Entity of a grid
// * @param[in] x Position in local coordinates with respect the entity
// * @param y Saturated conductivity at position x
// */
// void evaluate ( const typename Traits::ElementType& e,
// const typename Traits::DomainType& x,
// typename Traits::RangeType& y) const
// {
// typename GFWaterContent::Traits::RangeType water_content;
// typename GFWaterFlux::Traits::RangeType water_flux;
// // evaluate water flux
// _gf_water_flux->evaluate(e,x,water_flux);
// // evaluate water flux
// _gf_water_content->evaluate(e,x,water_content);
// // bind to entity
// _param->bind(e);
// // evaluate saturated conductivity
// auto y_raw = _param->effective_hydromechanic_dispersion_f()(water_flux,water_content);
// if constexpr (tensor_glyph)
// y = tensor_to_glyph(y_raw);
// else
// y = y_raw;
// }
// /*-------------------------------------------------------------------*//**
// * @brief Function that returns a grid view valid for this grid
// * function
// *
// * @return Grid view
// */
// const typename Traits::GridViewType& getGridView() const
// {
// return _gv;
// }
// private:
// const typename Traits::GridViewType _gv;
// const std::shared_ptr<const Parameter> _param;
// const std::shared_ptr<const GFWaterFlux> _gf_water_flux;
// const std::shared_ptr<const GFWaterContent> _gf_water_content;
// };
template<class T, class Parameter, class GFWaterFlux, class GFWaterContent>
template<class T, class Parameter, class GFWaterFlux, class GFWaterContent, int dimRange = 3>
class EffectiveHydrodynamicDispersionAdapter
: public Dune::PDELab::GridFunctionBase<
Dune::PDELab::GridFunctionTraits<
typename T::GridView,
typename T::RangeField,
T::dim,
Dune::FieldVector<typename T::RangeField,T::dim>
dimRange*dimRange,
Dune::FieldVector<typename T::RangeField,dimRange*dimRange>
>
,EffectiveHydrodynamicDispersionAdapter<T,Parameter,GFWaterFlux,GFWaterContent>
,EffectiveHydrodynamicDispersionAdapter<T,Parameter,GFWaterFlux,GFWaterContent,dimRange>
>
{
public:
using Traits = Dune::PDELab::GridFunctionTraits<
typename T::GridView,
typename T::RangeField,
T::dim,
Dune::FieldVector<typename T::RangeField,T::dim>
dimRange*dimRange,
Dune::FieldVector<typename T::RangeField,dimRange*dimRange>
>;
private:
......@@ -147,7 +57,7 @@ public:
* @param[in] x Position in local coordinates with respect the entity
* @param y Hydrodynamic dispersion at position x
*/
void inline evaluate ( const typename Traits::ElementType& e,
void evaluate ( const typename Traits::ElementType& e,
const typename Traits::DomainType& x,
typename Traits::RangeType& y) const
{
......@@ -161,10 +71,22 @@ public:
// bind to entity
_param->bind(e);
// evaluate hydrodynamic dispersion
auto y_raw = _param->hydrodynamic_dispersion_f()(water_flux,water_content);
auto y_tensor = _param->hydrodynamic_dispersion_f()(water_flux,water_content);
using Range = typename Traits::RangeType;
using TensorRange = decltype(y_tensor);
int count = 0;
if constexpr (dimRange!=TensorRange::rows or dimRange!=TensorRange::cols)
y = Range(0.);
for (int i = 0; i < T::dim; ++i)
y[i] = y_raw[i][i];
std::cout << y_tensor << std::endl;
for (int i = 0; i < dimRange; ++i)
for (int j = 0; j < dimRange; ++j, ++count)
if (i<TensorRange::rows and j<TensorRange::rows)
y[count] = y_tensor[i][j];
std::cout << y << std::endl;
}
/*-------------------------------------------------------------------*//**
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment