fract_thf_fv1_impl.h

Go to the documentation of this file.

/* 
 * SPDX-FileCopyrightText: 2025 Gesellschaft fuer Anlagen- und Reaktorsicherheit gGmbH
 * SPDX-License-Identifier: EUPL-1.2
 * SPDX-FileContributor: Dmitry Logashenko, Sabine Stichel
 * SPDX-FileContributor: Goethe Universität Frankfurt
 * SPDX-FileType: SOURCE
 * 
 * This file is part of d3f++.
 * d3f++ is an extension for UG4. Licensing information and citation requirements of UG4 are provided in LICENSES/UG4-LGPL_2.1
 */
 
 
/*
 * fract_thf_fv1_impl.h - class member implementation of the discretization of
 * the FV-discretization of the thermohaline flow in fractured porous media.
 *
 * Based on the modules by A. Vogel.
 */
 
/* UG4 headers: */
#include "lib_disc/spatial_disc/disc_util/geom_provider.h"
#include "lib_disc/spatial_disc/disc_util/fv1_geom.h"
#include "lib_disc/spatial_disc/user_data/const_user_data.h"
#include "lib_disc/spatial_disc/disc_util/conv_shape.h"
#ifdef UG_FOR_LUA
#include "bindings/lua/lua_user_data.h"
#endif
 
namespace ug{
namespace d3f{
 
//  Local discretization interface
 
template<typename TDomain>
void FractTHF_FV1<TDomain>::prepare_setting
(
  const std::vector<LFEID> & vLfeID,
  bool bNonRegular
)
{
//  check the grid
  if (bNonRegular)
    UG_THROW ("ERROR in FractTHF_FV1:"
        " The discretization does not support hanging nodes.\n");
 
//  check number of the components
  if (vLfeID.size () != 3)
    UG_THROW ("FractTHF_FV1: The density driven flow requires two coponents in every vertex.");
    /* one component for the real part and one for the imaginary part */
 
//  check whether these are the LagrangeP1 elements
  if (vLfeID[0] != LFEID(LFEID::LAGRANGE, dim, 1) || vLfeID[1] != LFEID(LFEID::LAGRANGE, dim, 1) || vLfeID[2] != LFEID(LFEID::LAGRANGE, dim, 1))
    UG_THROW ("FractTHF_FV1: This discretization works with the LagrangeP1-elements only.");
}
 
// assembling
 
template<typename TDomain>
void FractTHF_FV1<TDomain>::prep_assemble_loop ()
{
  if (m_spVolStabData.valid())
  //  set the "volume extensions" to 0
    m_spVolStabData->reset_flux();
 
  if (m_sss_mngr.valid ())
  {
  //  reset the markers of the singular sources and sinks
    m_sss_mngr->init_all_point_sss ();
    m_sss_mngr->init_all_line_sss ();
  }
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::prepare_element_loop
(
  ReferenceObjectID roid, 
  int si 
)
{
//  check whether we are in a degenerated fracture
  if (! m_spFractManager.valid())
    m_isFracture = false;
  else
  {
    if (! m_spFractManager->is_closed ())
      UG_THROW ("FractTHF_FV1: Fracture manager not closed");
    m_isFracture = m_spFractManager->contains (si);
  }
  
//  check the imports
  if (!m_imDensityCo.data_given())
    UG_THROW ("FractTHF_FV1: Missing Import 'density at corners'.");
  if (!m_imDensityIP.data_given())
    UG_THROW ("FractTHF_FV1: Missing Import 'density at IPs'.");
  if (!m_imViscosityCo.data_given())
    UG_THROW ("FractTHF_FV1: Missing Import 'viscosity at corners'.");
  if (!m_imViscosityIP.data_given())
    UG_THROW ("FractTHF_FV1: Missing Import 'viscosity at IP's'.");
  if (!m_imConstGravity.data_given())
    UG_THROW ("FractTHF_FV1: Missing Import 'gravity'.");
  if (!m_imPorosity.data_given())
    UG_THROW ("FractTHF_FV1: Missing Import 'porosity'.");
  if(m_spVolStabData.valid() && !m_imOldDensityCo.data_given())
    UG_THROW ("FractTHF_FV1: Missing Import 'density at corners'.");
  
//  get the gravity
  if (m_imConstGravity.constant())
    (*m_imConstGravity.user_data()) (m_Gravity, MathVector<dim>(), 0.0, 0);
  else
    UG_THROW ("FractTHF_FV1: Gravity must be constant.");
  
//  call the specific routine:
  if (! m_isFracture)
    this->template bulk_prepare_element_loop<TElem> (roid, si);
  else
    this->template fract_prepare_element_loop<TElem> (roid, si);
 
//  check if the volume stabilization is active
  m_bVolStabDataActive = (m_spVolStabData.valid() && m_spVolStabData->flux_is_active());
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::bulk_prepare_element_loop
(
  ReferenceObjectID roid, 
  int si 
)
{
  typedef FV1Geometry<TElem, dim> TBulkFVGeom;
  
//  check the imports
  if (!m_imPermeability.data_given())
    UG_THROW ("FractTHF_FV1: Missing Import 'full-dim. permeability'.");
  if (!m_imDiffusion.data_given())
    UG_THROW ("FractTHF_FV1: Missing Import 'full_dim. diffusion'.");
  if (!m_imThermalConductivity.data_given())
    UG_THROW ("FractTHF_FV1: Missing Import 'full_dim. thermal conductivity'.");
  
//  check, that upwind has been set
  if (m_spUpwind.invalid())
    UG_THROW("FractTHF_FV1: Upwind for the transport has not been set.");
  if (m_spUpwindEnergy.invalid())
    UG_THROW("FractTHF_FV1: Upwind for the energy has not been set.");
 
  TBulkFVGeom& bulk_geo = GeomProvider<TBulkFVGeom>::get();
  static const int refDim = TElem::dim;
  
//  set local IP coordinates
  const MathVector<refDim>* vSCVFip = bulk_geo.scvf_local_ips();
  size_t numSCVFip = bulk_geo.num_scvf_ips();
  m_imDensityIP.template set_local_ips<refDim> (vSCVFip, numSCVFip, false);
  m_imViscosityIP.template set_local_ips<refDim> (vSCVFip, numSCVFip, false);
  
  m_imFractDarcyVelIP.template set_local_ips<refDim> (vSCVFip, numSCVFip, false);//to delete
 
  const MathVector<refDim>* vSCVip = bulk_geo.scv_local_ips();
  size_t numSCVip = bulk_geo.num_scv_ips();
  m_imDensityCo.template set_local_ips<refDim> (vSCVip, numSCVip, false);
  m_imViscosityCo.template set_local_ips<refDim> (vSCVip, numSCVip, false);
  if (m_spVolStabData.valid ())
    m_imOldDensityCo.template set_local_ips<refDim> (vSCVip, numSCVip, 1, false);
  
  const MathVector<refDim>* coe_local = bulk_geo.coe_local ();
  m_imPorosity.template set_local_ips<refDim> (coe_local, 1, false);
  m_imPermeability.template set_local_ips<refDim> (coe_local, 1, false);
  m_imDiffusion.template set_local_ips<refDim> (coe_local, 1, false);
  m_imThermalConductivity.template set_local_ips<refDim> (coe_local, 1, false);
  
//  init upwind for element type
  if(! m_spUpwind->template set_geometry_type<TBulkFVGeom>(bulk_geo))
    UG_THROW("FractTHF_FV1: Cannot init upwind for the transport for bulk element type.");
//  init upwind for element type
  if(! m_spUpwindEnergy->template set_geometry_type<TBulkFVGeom>(bulk_geo))
    UG_THROW("FractTHF_FV1: Cannot init upwind for the energy for bulk element type.");
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::fract_prepare_element_loop
(
  ReferenceObjectID roid, 
  int si 
)
{
  typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
  static const int refDim = TElem::dim; // dimensionality of the element, not the side!
  ref_elem_type& rRefElem = Provider<ref_elem_type>::get ();
  
//  initialize the pointer to the FV geometry for fracture elements
  m_pFractGeo = &GeomProvider<TFractFVGeom>::get (LFEID (LFEID::LAGRANGE, low_dim, 1), 1);
  
//  check the imports
  if (!m_imAperture.data_given())
    UG_THROW ("FractTHF_FV1: Missing Import 'fracture width (aperture)'.");
  if (!m_imFractPermeability.data_given())
    UG_THROW ("FractTHF_FV1: Missing Import 'permeability along fracture'.");
  if (!m_imOrthoPermeability.data_given())
    UG_THROW ("FractTHF_FV1: Missing Import 'permeability of the bulk-fracture interface'.");
  if (!m_imFractDiffusion.data_given())
    UG_THROW ("FractTHF_FV1: Missing Import 'diffusion along fracture'.");
  if (!m_imOrthoDiffusion.data_given())
    UG_THROW ("FractTHF_FV1: Missing Import 'diffusion through the bulk-fracture interface'.");
  if (!m_imFractThermCond.data_given())
    UG_THROW ("FractTHF_FV1: Missing Import 'thermal conductivity along fracture'.");
  if (!m_imOrthoThermCond.data_given())
    UG_THROW ("FractTHF_FV1: Missing Import 'thermal conductivity through the bulk-fracture interface'.");
  if(!m_imFractDarcyVelIP.data_given())
    UG_THROW("FractTHF_FV1: Missing Import: Darcy Velocity in fracture.");
  
//  check, that upwind has been set
  if (m_spUpwind.invalid())
    UG_THROW("FractTHF_FV1: Upwind for the transport has not been set.");
//  check, that upwind has been set
  if (m_spUpwindEnergy.invalid())
    UG_THROW("FractTHF_FV1: Upwind for the energy has not been set.");
  
//  set up local ip coordinates for corner import parameters
//  REMARK: Note that for the fracture elements, values of the corner import
//  parameters are indexed not by scv (as for the normal elements) but by
//  the corner indices in the reference element
  m_imDensityCo.template set_local_ips<refDim> (rRefElem.corners(), ref_elem_type::numCorners, false);
  m_imViscosityCo.template set_local_ips<refDim> (rRefElem.corners(), ref_elem_type::numCorners, false);
  if(m_spVolStabData.valid())
    m_imOldDensityCo.template set_local_ips<refDim> (rRefElem.corners(), ref_elem_type::numCorners, 1, false);
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::finish_element_loop ()
{
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::prepare_element
(
  const LocalVector & u, 
  GridObject * elem, 
  ReferenceObjectID roid,  // id of reference element used for assembling
  const position_type vCornerCoords [] 
)
{
  TElem * pElem = static_cast<TElem*> (elem);
  
//  call the specific routine:
  if (! m_isFracture)
    this->template bulk_prepare_element<TElem> (u, pElem, vCornerCoords);
  else
    this->template fract_prepare_element<TElem> (u, pElem, vCornerCoords);
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::bulk_prepare_element
(
  const LocalVector & u, 
  TElem * elem, 
  const position_type vCornerCoords [] 
)
{
  typedef FV1Geometry<TElem, dim> TBulkFVGeom;
  
//  update the FV geometry for this element
  TBulkFVGeom& bulk_geo = GeomProvider<TBulkFVGeom>::get();
  try
  {
    bulk_geo.update (elem, vCornerCoords, &(this->subset_handler()));
  }
  UG_CATCH_THROW("FractTHF_FV1: Cannot update the Finite Volume Geometry for a bulk element.");
 
//  set global positions for user data
  const MathVector<dim>* vSCVFip = bulk_geo.scvf_global_ips();
  size_t numSCVFip = bulk_geo.num_scvf_ips();
  m_imDensityIP.set_global_ips (vSCVFip, numSCVFip);
  m_imViscosityIP.set_global_ips (vSCVFip, numSCVFip);
  
  m_imFractDarcyVelIP.set_global_ips (vSCVFip, numSCVFip);//to delete
 
  const MathVector<dim>* vSCVip = bulk_geo.scv_global_ips();
  size_t numSCVip = bulk_geo.num_scv_ips();
  m_imDensityCo.set_global_ips (vSCVip, numSCVip);
  m_imViscosityCo.set_global_ips (vSCVip, numSCVip);
  if(m_spVolStabData.valid())
    m_imOldDensityCo.set_global_ips (vSCVip, numSCVip);
  
  const MathVector<dim>* coe_global = bulk_geo.coe_global ();
  m_imPorosity.set_global_ips (coe_global, 1);
  m_imPermeability.set_global_ips (coe_global, 1);
  m_imDiffusion.set_global_ips (coe_global, 1);
  m_imThermalConductivity.set_global_ips (coe_global, 1);
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::fract_prepare_element
(
  const LocalVector & u, 
  TElem * elem, 
  const position_type vCornerCoords [] 
)
{
  typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
  ref_elem_type& rRefElem = Provider<ref_elem_type>::get ();
  
//  get the non-degenerated sides of the fracture element
  try
  {
    m_spFractManager->get_layer_sides
      (elem,
        m_numFractCo, m_innerFractSide, m_innerFractSideIdx, m_innerSideCo,
        m_outerFractSide, m_outerFractSideIdx, m_outerSideCo,
        m_assCo);
  }
  UG_CATCH_THROW("FractTHF_FV1: Cannot find orientation of a fracture element.");
  
//  compute the FV geometry of the inner side
  MathVector<dim> vSideCornerCoords [maxFractSideCorners];
  try
  {
    for (size_t co = 0; co < m_numFractCo; co++)
      vSideCornerCoords [co] = vCornerCoords [m_innerSideCo [co]];
    m_pFractGeo->update (m_innerFractSide, vSideCornerCoords, &(this->subset_handler()));
  }
  UG_CATCH_THROW("FractTHF_FV1: Cannot update the Finite Volume Geometry for a fracture element.");
  size_t numSCVFip = m_pFractGeo->num_scvf_ips ();
  
//  convert local coordinates of the side into the local coordinates of the element (for the input parameters)
  MathVector<dim> vSideLocCornerCoords [maxFractSideCorners];
  try
  {
    DimReferenceMapping<low_dim, dim>& rMapping
      = ReferenceMappingProvider::get<low_dim, dim> (m_innerFractSide->reference_object_id ());
    for (size_t co = 0; co < m_numFractCo; co++)
      vSideLocCornerCoords [co] = rRefElem.corner (m_innerSideCo [co]);
    rMapping.update (vSideLocCornerCoords);
    
    rMapping.local_to_global (m_elem_loc_coe, *(m_pFractGeo->coe_local ()));
    rMapping.local_to_global (m_elem_loc_scvf, m_pFractGeo->scvf_local_ips (), numSCVFip);
  }
  UG_CATCH_THROW("FractTHF_FV1: Cannot transform local side coordinates to local element coordinates in a fracture element.");
  
//  set local IP coordinates
  m_imDensityIP.template set_local_ips<dim> (m_elem_loc_scvf, numSCVFip);
  m_imViscosityIP.template set_local_ips<dim> (m_elem_loc_scvf, numSCVFip);
  
  m_imAperture.template set_local_ips<dim> (&m_elem_loc_coe, 1);
  m_imPorosity.template set_local_ips<dim> (&m_elem_loc_coe, 1);
  m_imFractPermeability.template set_local_ips<dim> (&m_elem_loc_coe, 1);
  m_imOrthoPermeability.template set_local_ips<dim> (&m_elem_loc_coe, 1);
  m_imFractDiffusion.template set_local_ips<dim> (&m_elem_loc_coe, 1);
  m_imOrthoDiffusion.template set_local_ips<dim> (&m_elem_loc_coe, 1);
  m_imFractThermCond.template set_local_ips<dim> (&m_elem_loc_coe, 1);
  m_imOrthoThermCond.template set_local_ips<dim> (&m_elem_loc_coe, 1);
  
//  set global positions for user data
  const MathVector<dim>* vSCVFip = m_pFractGeo->scvf_global_ips();
  m_imDensityIP.set_global_ips (vSCVFip, numSCVFip);
  m_imViscosityIP.set_global_ips (vSCVFip, numSCVFip);
  
  const MathVector<dim>* coe_global = m_pFractGeo->coe_global ();
  m_imAperture.set_global_ips (coe_global, 1);
  m_imPorosity.set_global_ips (coe_global, 1);
  m_imFractPermeability.set_global_ips (coe_global, 1);
  m_imOrthoPermeability.set_global_ips (coe_global, 1);
  m_imFractDiffusion.set_global_ips (coe_global, 1);
  m_imOrthoDiffusion.set_global_ips (coe_global, 1);
  m_imFractThermCond.set_global_ips (coe_global, 1);
  m_imOrthoThermCond.set_global_ips (coe_global, 1);
  
//  set up global ip coordinates for corner import parameters
//  REMARK: Note that for the fracture elements, values of the corner import
//  parameters are indexed not by scv (as for the normal elements) but by
//  the corner indices in the reference element
  m_imDensityCo.set_global_ips (vCornerCoords, ref_elem_type::numCorners);
  m_imViscosityCo.set_global_ips (vCornerCoords, ref_elem_type::numCorners);
  if(m_spVolStabData.valid())
    m_imOldDensityCo.set_global_ips (vCornerCoords, ref_elem_type::numCorners);
 
//  set global and local positions for velocity import
  m_imFractDarcyVelIP.template set_local_ips<dim> (m_elem_loc_scvf, numSCVFip);
  m_imFractDarcyVelIP.set_global_ips(vSCVFip, numSCVFip);
 
//  compute the orthogonal gravity
  MathVector<dim> outerNormal;
  number outerNormalNorm;
  SideNormal<ref_elem_type, dim> (outerNormal, m_outerFractSideIdx, vCornerCoords);
  if ((outerNormalNorm = VecLength (outerNormal)) < 1e-32)
    UG_THROW ("FractTHF_FV1: Cannot get the normal to a fracture.")
  m_orthGravity = VecDot (m_Gravity, outerNormal) / outerNormalNorm;
  
//  init upwind for element type
  if(! m_spUpwind->template set_geometry_type<TFractFVGeom>(*m_pFractGeo))
    UG_THROW("FractTHF_FV1: Cannot init upwind for the transport for fracture element type.");
//  init upwind for element type
  if(! m_spUpwindEnergy->template set_geometry_type<TFractFVGeom>(*m_pFractGeo))
    UG_THROW("FractTHF_FV1: Cannot init upwind for the energy for fracture element type.");
}
 
template<typename TDomain>
template<typename TFVGeom, typename TConsGravity>
void FractTHF_FV1<TDomain>::compute_ip_Darcy_velocity
(
  MathVector<dim>& Vel, 
  size_t ip, 
  const TFVGeom& geo, 
    TConsGravity& ConsGravityMethod, 
    MathVector<TFVGeom::dim> vConsGravity[], 
    number vPressure[], 
  number Viscosity 
)
{
  const typename TFVGeom::SCVF& scvf = geo.scvf(ip);
 
//  Compute rho * g (as the consistent gravity force)
  ConsGravityMethod.template compute<dim>
    (Vel, scvf.local_ip(), scvf.JTInv(), scvf.local_grad_vector(), vConsGravity);
 
// The pressure-gradient part:
  for (size_t sh = 0; sh < scvf.num_sh(); sh++)
    VecScaleAppend (Vel, - vPressure[sh], scvf.global_grad(sh));
 
// The viscosity factor:
  VecScale (Vel, Vel, 1 / Viscosity);
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::ass_dA_elem
(
  LocalVector & d, 
  const LocalVector & u, 
  GridObject * elem, 
  const position_type vCornerCoords [] 
)
{
  TElem * pElem = static_cast<TElem*> (elem);
 
//  call the specific routine:
  if (! m_isFracture)
  {
    this->template bulk_ass_dA_elem<TElem> (d, u, pElem, vCornerCoords);
  }
  else
  {
    this->template fract_ass_dA_elem<TElem> (d, u, pElem, vCornerCoords);
    this->template fract_bulk_ass_dA_elem<TElem> (d, u, pElem, vCornerCoords);
  }
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::bulk_ass_dA_elem
(
  LocalVector & d, 
  const LocalVector & u, 
  TElem * pElem, 
  const position_type vCornerCoords [] 
)
{
  typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
  typedef FV1Geometry<TElem, dim> TBulkFVGeom;
  typedef StdLinConsistentGravityX<ref_elem_type::dim> TConsGravity;
 
  MathMatrix<dim,dim> Diffusion[TBulkFVGeom::numSCVF];
  MathMatrix<dim,dim> ThermCond[TBulkFVGeom::numSCVF];
  MathVector<dim> Dgrad_c_ip, Dgrad_T_ip;
  MathVector<dim> grad_c_ip, grad_T_ip;
  MathVector<dim> Vel[TBulkFVGeom::numSCVF], DarcyVel[TBulkFVGeom::numSCVF];
 
//  Get finite volume geometry
  const TBulkFVGeom& bulk_geo = GeomProvider<TBulkFVGeom>::get();
 
    const size_t numSh = bulk_geo.num_sh();
  const size_t numScvf = bulk_geo.num_scvf();
 
//  Consistent gravity and its derivative at corners
  TConsGravity ConsGravityMethod;
    MathVector<ref_elem_type::dim> vConsGravity [ref_elem_type::numCorners];
 
//  Prepare the consistent gravity
  try
  {
    ConsGravityMethod.template prepare<dim>
      (vConsGravity, numSh, vCornerCoords, m_imDensityCo.values(), m_Gravity);
  }
  UG_CATCH_THROW ("FractTHF_FV1::bulk_ass_dA_elem: Cannot prepare Consistent Gravity.");
  
//  Get the corner pressure
  number vPressure [ref_elem_type::numCorners];
  for (size_t sh = 0; sh < numSh; sh++)
    vPressure[sh] = u(_P_, sh);
  
//  Compute Diffusion - Dispersion
  MatScale(Diffusion[0], m_imPorosity[0], m_imDiffusion[0]);
  for (size_t ip = 1; ip < numScvf; ip++)
    Diffusion[ip] = Diffusion[0];
  //TODO add Dispersion
  
//  Compute the thermal conductivity
  ThermCond[0] = m_imThermalConductivity[0];
  for (size_t ip = 1; ip < numScvf; ip++)
    ThermCond[ip] = ThermCond[0];
  
//  Compute the Darcy velocity at the ips and the convection shapes
  for (size_t ip = 0; ip < numScvf; ip++)
  {
      this->template compute_ip_Darcy_velocity<TBulkFVGeom, TConsGravity>
          (Vel[ip], ip, bulk_geo, ConsGravityMethod, vConsGravity, vPressure,
            m_imViscosityIP[ip]);
      MatVecMult (DarcyVel[ip], m_imPermeability[0], Vel[ip]);
  }
  
//  compute upwind shapes for transport equation
  if(!m_spUpwind->update(&bulk_geo, DarcyVel, Diffusion, false))
    UG_THROW("FractTHF_FV1::bulk_ass_dA_elem: Cannot compute convection shapes for the transport equation.");
 
//  compute upwind shapes for transport equation
  if(!m_spUpwindEnergy->update(&bulk_geo, DarcyVel, ThermCond, false))
    UG_THROW("FractTHF_FV1::bulk_ass_dA_elem: Cannot compute convection shapes for the energy equation.");
 
//  get a const (!!) reference to the upwind
  const IConvectionShapes<dim>& convShape
    = *const_cast<const IConvectionShapes<dim>*>(m_spUpwind.get());
  
//  get a const (!!) reference to the upwind
  const IConvectionShapes<dim>& convShapeT
    = *const_cast<const IConvectionShapes<dim>*>(m_spUpwindEnergy.get());
  
  for (size_t ip = 0; ip < numScvf; ip++)
  {
  //  Get current SCVF
    const typename TBulkFVGeom::SCVF& scvf = bulk_geo.scvf(ip);
 
    // The the gradient of c and T at the ip:
    VecSet (grad_c_ip, 0.0); VecSet (grad_T_ip, 0.0);
    for (size_t sh = 0; sh < scvf.num_sh(); sh++)
    {
      VecScaleAppend (grad_c_ip, u(_C_, sh), scvf.global_grad(sh));
      VecScaleAppend (grad_T_ip, u(_T_, sh), scvf.global_grad(sh));
    }
 
  // Transport Equation
 
    number flux = 0;
    
  //  Compute Convective Flux
    for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
      flux += convShape(ip, sh) * u(_C_, sh);
 
  //  Compute Diffusive Flux
    MatVecMult(Dgrad_c_ip, Diffusion[ip], grad_c_ip);
    const number diffFlux = VecDot(Dgrad_c_ip, scvf.normal());
 
  //  Sum total flux
    flux -= diffFlux;
    if(!m_BoussinesqTransport) flux *= m_imDensityIP[ip];
 
  //  Add contribution to transport equation
    d(_C_,scvf.from()) += flux;
    d(_C_,scvf.to()) -= flux;
 
  // Flow Equation
 
  //  Compute flux
    flux = VecDot(DarcyVel[ip], scvf.normal());
    if(!m_BoussinesqFlow) flux *= m_imDensityIP[ip];
 
  //  Add contribution to flow equation
    d(_P_,scvf.from()) += flux;
    d(_P_,scvf.to()) -= flux;
    
  //  Volume stabilization: Add the continuity equation to the transport equation
    if(m_bVolStabDataActive)
    {
      m_spVolStabData->stiff(pElem->vertex(scvf.from())) -= flux;
      m_spVolStabData->stiff(pElem->vertex(scvf.to())) += flux;
    }
 
  // Energy Equation
 
  //  Compute Convective Flux
    flux = 0;
    for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
      flux += convShapeT(ip, sh) * u(_T_, sh);
    flux *= m_imHeatCapacityFluid;
 
    if(m_BoussinesqEnergy)  flux *= m_BoussinesqDensity;
    else          flux *= m_imDensityIP[ip];
 
  //  Compute Diffusive Flux
    MatVecMult(Dgrad_T_ip, ThermCond[ip], grad_T_ip);
    const number thermDiffFlux = VecDot(Dgrad_T_ip, scvf.normal());
 
  //  Sum total flux
    flux -= thermDiffFlux;
 
  //  Add contribution to transport equation
    d(_T_,scvf.from()) += flux;
    d(_T_,scvf.to()) -= flux;
  }
  
//  Assemble the singular sources and sinks
    if (m_sss_mngr.valid () && (m_sss_mngr->num_points () != 0 || m_sss_mngr->num_lines () != 0))
    {
      typedef typename domain_type::position_accessor_type t_pos_accessor;
      typedef typename THSingularSourcesAndSinks<dim>::template
        point_iterator<TElem,t_pos_accessor,TBulkFVGeom> t_pnt_sss_iter;
      typedef typename THSingularSourcesAndSinks<dim>::template
        line_iterator<TElem,t_pos_accessor,TBulkFVGeom> t_lin_sss_iter;
      
    t_pos_accessor& aaPos = this->domain()->position_accessor ();
    Grid& grid = (Grid&) *this->domain()->grid ();
    
    for(size_t ip = 0; ip < bulk_geo.num_scv(); ip++)
    {
      size_t co = bulk_geo.scv(ip).node_id ();
      
    //  point sources
      for (t_pnt_sss_iter pnt (m_sss_mngr.get (), pElem, grid, aaPos, bulk_geo, co);
        ! pnt.is_over (); ++pnt)
      {
        FVPointSourceOrSink<dim, th_point_sss_data<dim> > * pnt_sss = *pnt;
        if (! pnt_sss->marked_for (pElem, co))
          continue;
        pnt_sss->compute (pnt_sss->position (), this->time (), -1); //TODO: set the subset id instead of -1
        ass_sss_dA_elem (d, u, pElem, co, pnt_sss->intensity (), pnt_sss->concentration (), pnt_sss->temperature ());
      }
      
    //  line sources
      for (t_lin_sss_iter line (m_sss_mngr.get (), pElem, grid, aaPos, bulk_geo, co);
        ! line.is_over (); ++line)
      {
        FVLineSourceOrSink<dim, th_line_sss_data<dim> > * line_sss = *line;
        number len = VecDistance (line.seg_start (), line.seg_end ());
        line_sss->compute (line.seg_start (), this->time (), -1); //TODO: set the subset id instead of -1
        ass_sss_dA_elem (d, u, pElem, co, line_sss->intensity (), line_sss->concentration (), line_sss->temperature (), len);
      }
    }
  }
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::fract_ass_dA_elem
(
  LocalVector & d, 
  const LocalVector & u, 
  TElem * pElem, 
  const position_type vCornerCoords [] 
)
{
  const number half_fr_width = m_imAperture[0] / 2;
  
  MathMatrix<dim,dim> Diffusion[TFractFVGeom::maxNumSCVF];
  MathMatrix<dim,dim> ThermCond[TFractFVGeom::maxNumSCVF];
  MathVector<dim> Dgrad_c_ip, Dgrad_T_ip;
  MathVector<dim> grad_c_ip, grad_T_ip;
 
    const size_t numSh = m_pFractGeo->num_sh();
  const size_t numScvf = m_pFractGeo->num_scvf();
 
//  Get the corner values
  number vConcentration [maxFractSideCorners];
  number vTemperature [maxFractSideCorners];
  for (size_t sh = 0; sh < numSh; sh++)
  {
    size_t co = m_innerSideCo[sh];
    vConcentration[sh] = u(_C_, co);
    vTemperature[sh] = u(_T_, co);
  }
  
//  Compute Diffusion - Dispersion
  MatSet (Diffusion[0], 0);
  MatDiagSet (Diffusion[0], m_imPorosity[0] * m_imFractDiffusion[0]);
  for (size_t ip = 1; ip < numScvf; ip++)
    Diffusion[ip] = Diffusion[0];
  //TODO add Dispersion
  
//  Compute the thermal conductivity
  MatSet (ThermCond[0], 0);
  MatDiagSet (ThermCond[0], m_imFractThermCond[0]);
  for (size_t ip = 1; ip < numScvf; ip++)
    ThermCond[ip] = ThermCond[0];
  
//  compute upwind shapes for transport equation
  if(!m_spUpwind->update(m_pFractGeo, m_imFractDarcyVelIP.values(), Diffusion, false))
    UG_THROW("FractTHF_FV1::fract_ass_dA_elem: Cannot compute convection shapes for the transport equation.");
 
  const IConvectionShapes<dim>& convShape
    = *const_cast<const IConvectionShapes<dim>*>(m_spUpwind.get());
  
//  compute upwind shapes for energy equation
  if(!m_spUpwindEnergy->update(m_pFractGeo, m_imFractDarcyVelIP.values(), ThermCond, false))
    UG_THROW("FractTHF_FV1::fract_ass_dA_elem: Cannot compute convection shapes for the energy equation.");
 
  const IConvectionShapes<dim>& convShapeT
    = *const_cast<const IConvectionShapes<dim>*>(m_spUpwindEnergy.get());
  
//  Loop Sub Control Volume Faces (SCVF)
  for (size_t ip = 0; ip < numScvf; ip++)
  {
  //  Get current SCVF
    const typename TFractFVGeom::SCVF& scvf = m_pFractGeo->scvf(ip);
 
    // The ip-concentration and the gradient of c at the ip:
    VecSet (grad_c_ip, 0.0); VecSet (grad_T_ip, 0.0);
    for (size_t sh = 0; sh < scvf.num_sh(); sh++)
    {
      VecScaleAppend (grad_c_ip, vConcentration[sh], scvf.global_grad(sh));
      VecScaleAppend (grad_T_ip, vTemperature[sh], scvf.global_grad(sh));
    }
 
  // Transport Equation
 
    number flux = 0;
  
  //  Compute Convective Flux
    for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
      flux += convShape(ip, sh) * vConcentration[sh];
 
  //  Compute Diffusive Flux
    MatVecMult(Dgrad_c_ip, Diffusion[ip], grad_c_ip);
    const number diffFlux = VecDot(Dgrad_c_ip, scvf.normal());
 
  //  Sum total flux
    flux = (flux - diffFlux) * half_fr_width;
    if(!m_BoussinesqTransport) flux *= m_imDensityIP[ip];
 
  //  Add contribution to transport equation
    d(_C_, m_innerSideCo[scvf.from()]) += flux;
    d(_C_, m_innerSideCo[scvf.to()]) -= flux;
 
  // Flow Equation
 
  //  Compute flux
    flux = VecDot(m_imFractDarcyVelIP[ip], scvf.normal()) * half_fr_width;
 
    if(!m_BoussinesqFlow) flux *= m_imDensityIP[ip];
 
  //  Add contribution to flow equation
    d(_P_, m_innerSideCo[scvf.from()]) += flux;
    d(_P_, m_innerSideCo[scvf.to()]) -= flux;
    
  //  Volume stabilization: Add the continuity equation to the transport equation
    if(m_bVolStabDataActive)
    {
      m_spVolStabData->stiff(pElem->vertex(m_innerSideCo[scvf.from()])) -= flux;
      m_spVolStabData->stiff(pElem->vertex(m_innerSideCo[scvf.to()])) += flux;
    }
 
  // Energy Equation
 
  //  Compute Convective Flux
    flux = 0;
    for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
      flux += convShapeT(ip, sh) * vTemperature[sh];
    flux *= m_imHeatCapacityFluid;
 
    if(m_BoussinesqEnergy)  flux *= m_BoussinesqDensity;
    else          flux *= m_imDensityIP[ip];
 
  //  Compute Diffusive Flux
    MatVecMult(Dgrad_T_ip, ThermCond[ip], grad_T_ip);
    const number thermDiffFlux = VecDot(Dgrad_T_ip, scvf.normal());
 
  //  Sum total flux
    flux = (flux - thermDiffFlux) * half_fr_width;
 
  //  Add contribution to transport equation
    d(_T_, m_innerSideCo[scvf.from()]) += flux;
    d(_T_, m_innerSideCo[scvf.to()]) -= flux;
  }
  
//  Assemble the singular sources and sinks
    if (m_sss_mngr.valid () && m_sss_mngr->num_lines () != 0)
    {
      typedef typename domain_type::position_accessor_type t_pos_accessor;
      typedef typename THSingularSourcesAndSinks<dim>::template
        line_iterator<side_type,t_pos_accessor,TFractFVGeom> t_lin_sss_iter;
      
    t_pos_accessor& aaPos = this->domain()->position_accessor ();
    Grid& grid = (Grid&) *this->domain()->grid ();
    
    for(size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
    {
    //  Get the corner of the face
      size_t side_co = m_pFractGeo->scv(ip).node_id ();
    //  Get associated node of the element (not side!)
      size_t co = m_innerSideCo [m_pFractGeo->scv(ip).node_id ()];
      
    //  line sources (that correspond to the point sources)
      for (t_lin_sss_iter line (m_sss_mngr.get (), m_innerFractSide, grid, aaPos, *m_pFractGeo, side_co);
        ! line.is_over (); ++line)
      {
        FVLineSourceOrSink<dim, th_line_sss_data<dim> > * line_sss = *line;
        if (! line_sss->marked_for (m_innerFractSide, side_co))
          continue;
        line_sss->compute (line.seg_start (), this->time (), -1); //TODO: set the subset id instead of -1
        ass_sss_dA_elem (d, u, pElem, co, line_sss->intensity (), line_sss->concentration (), line_sss->temperature (), 0.5);
        /* Remark: "0.5" because the source is taken into account twise. */
      }
    }
  }
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::fract_bulk_ass_dA_elem
(
  LocalVector & d, 
  const LocalVector & u, 
  TElem * pElem, 
  const position_type vCornerCoords [] 
)
{
  const number orthDiffusion = m_imPorosity[0] * m_imOrthoDiffusion[0];
  
  const number half_fr_width = m_imAperture[0] / 2;
  
//  loop over the corners of the inner side
  for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
  {
  //  Get current SCV
    const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv(ip);
 
  //  Get associated node of the element (not side!)
    const int co = m_innerSideCo [scv.node_id()];
    
  //  Get the corner values
    const number orthC_f = u(_C_, co);
    const number orthC_m = u(_C_, m_assCo[co]);
    const number orthP_f = u(_P_, co);
    const number orthP_m = u(_P_, m_assCo[co]);
    const number orthT_f = u(_T_, co);
    const number orthT_m = u(_T_, m_assCo[co]);
    const number fractDensity = m_imDensityCo[co];
    const number orthDensity = m_imDensityCo [m_assCo[co]];
    const number orthViscosity = m_imViscosityCo [m_assCo[co]];
  
  //  Flux velocity between the fracture and the bulk medium
    const number orthVelocity = (m_orthGravity * (orthDensity - fractDensity)
                    - (orthP_m - orthP_f) / half_fr_width)
                      * m_imOrthoPermeability[0] / orthViscosity;
 
    number flux;
 
  // Transport Equation
 
        flux = orthDiffusion * (orthC_m - orthC_f) / half_fr_width;
        /* We use the full upwind here: */
        flux -= orthVelocity * ((orthVelocity >= 0)? orthC_f : orthC_m);
        flux *= scv.volume();
        if (! m_BoussinesqTransport) flux *= orthDensity;
        d(_C_, m_assCo[co]) += flux;
        d(_C_, co) -= flux;
 
  // Flow Equation
 
        flux = orthVelocity * scv.volume();
        if (! m_BoussinesqFlow) flux *= orthDensity;
        d(_P_, m_assCo[co]) -= flux;
        d(_P_, co) += flux;
        
    //  Volume stabilization: Add the continuity equation to the transport equation
    if(m_bVolStabDataActive)
    {
      m_spVolStabData->stiff(pElem->vertex(m_assCo[co])) += flux;
      m_spVolStabData->stiff(pElem->vertex(co)) -= flux;
    }
 
  // Energy Equation
 
        /* We use the full upwind here: */
        flux = - orthVelocity * ((orthVelocity >= 0)? orthT_f : orthT_m)
            * m_imHeatCapacityFluid
            * (m_BoussinesqEnergy? m_BoussinesqDensity : orthDensity);
        flux += m_imOrthoThermCond[0] * (orthT_m - orthT_f) / half_fr_width;
        flux *= scv.volume();
        d(_T_, m_assCo[co]) += flux;
        d(_T_, co) -= flux;
  }
}
 
template<typename TDomain>
template<typename TFVGeom, typename TConsGravity, size_t maxCorners>
void FractTHF_FV1<TDomain>::compute_J_ip_Darcy_velocity
(
  MathVector<dim>& Vel, 
  MathVector<dim> Vel_c[], 
  MathVector<dim> Vel_p[], 
  MathVector<dim> Vel_T[], 
  size_t ip, 
  const TFVGeom& geo, 
    TConsGravity& ConsGravityMethod, 
    MathVector<TFVGeom::dim> vConsGravity_c[][maxCorners], 
    MathVector<TFVGeom::dim> vConsGravity_T[][maxCorners], 
    number vPressure[], 
  number Viscosity, 
  number Viscosity_c[], 
  number Viscosity_T[] 
)
{
  const typename TFVGeom::SCVF& scvf = geo.scvf(ip);
  const number InvVisco = 1 / Viscosity;
    const size_t numSh = geo.num_sh();
    
//  Compute the derivatives of rho g - grad p, multiplied by the viscosity
    for (size_t sh = 0; sh < numSh; sh++)
    {
      ConsGravityMethod.template compute<dim>
      (Vel_c[sh], scvf.local_ip(), scvf.JTInv(), scvf.local_grad_vector(),
        vConsGravity_c[sh]);
    VecScale (Vel_c[sh], Vel_c[sh], InvVisco);
    
      ConsGravityMethod.template compute<dim>
      (Vel_T[sh], scvf.local_ip(), scvf.JTInv(), scvf.local_grad_vector(),
        vConsGravity_T[sh]);
    VecScale (Vel_T[sh], Vel_T[sh], InvVisco);
    
    VecScale (Vel_p[sh], scvf.global_grad(sh), -InvVisco);
    }
    
//  Take into account the derivative of the viscosity wrt c
  if (Viscosity_c != NULL)
    for (size_t sh = 0; sh < numSh; sh++)
      VecScaleAppend (Vel_c[sh], -Viscosity_c[sh] * InvVisco, Vel);
    
//  Take into account the derivative of the viscosity wrt T
  if (Viscosity_T != NULL)
    for (size_t sh = 0; sh < numSh; sh++)
      VecScaleAppend (Vel_T[sh], -Viscosity_T[sh] * InvVisco, Vel);
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::ass_JA_elem
(
  LocalMatrix & J, 
  const LocalVector & u, 
  GridObject * elem, 
  const position_type vCornerCoords [] 
)
{
  TElem * pElem = static_cast<TElem*> (elem);
 
//  call the specific routine:
  if (! m_isFracture)
    this->template bulk_ass_JA_elem<TElem> (J, u, pElem, vCornerCoords);
  else
  {
    this->template fract_ass_JA_elem<TElem> (J, u, pElem, vCornerCoords);
    this->template fract_bulk_ass_JA_elem<TElem> (J, u, pElem, vCornerCoords);
  }
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::bulk_ass_JA_elem
(
  LocalMatrix & J, 
  const LocalVector & u, 
  TElem * pElem, 
  const position_type vCornerCoords [] 
)
{
  typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
  typedef FV1Geometry<TElem, dim> TBulkFVGeom;
  typedef StdLinConsistentGravityX<ref_elem_type::dim> TConsGravity;
  static const size_t numCorners = ref_elem_type::numCorners;
 
  //  Get finite volume geometry
  const TBulkFVGeom& bulk_geo = GeomProvider<TBulkFVGeom>::get();
 
  const size_t numSh = bulk_geo.num_sh();
  const size_t numScvf = bulk_geo.num_scvf();
 
  MathMatrix<dim,dim> Diffusion[TBulkFVGeom::numSCVF];
  MathMatrix<dim,dim> ThermCond[TBulkFVGeom::numSCVF];
  MathVector<dim> Dgrad;
  MathVector<dim> Vel[TBulkFVGeom::numSCVF], DarcyVel[TBulkFVGeom::numSCVF];
  MathVector<dim> Vel_c[numCorners], Vel_p[numCorners], Vel_T[numCorners];
  MathVector<dim> vDDarcyVel_c[numCorners], vDDarcyVel_p[numCorners], vDDarcyVel_T[numCorners];
  MathVector<dim> grad_c_ip, grad_T_ip;
  number vDFlux_c [numCorners], vDFlux_p [numCorners], vDFlux_T [numCorners];
  number vViscosity_c [numCorners], * pViscosity_c;
  number vViscosity_T [numCorners], * pViscosity_T;
 
//  Consistent gravity and its derivative at corners
  TConsGravity ConsGravityMethod;
    MathVector<ref_elem_type::dim> vConsGravity [numCorners];
    MathVector<ref_elem_type::dim> vConsGravity_c [numCorners][numCorners];
    MathVector<ref_elem_type::dim> vConsGravity_T [numCorners][numCorners];
 
//  Prepare the consistent gravity
  try
  {
    ConsGravityMethod.template prepare<dim>
      (vConsGravity, numSh, vCornerCoords, m_imDensityCo.values(), m_Gravity);
    
    number co_density [numCorners];
    memset (co_density, 0, ref_elem_type::numCorners * sizeof (number));
    for (size_t sh = 0; sh < numSh; sh++)
    { //TODO: use scv, not merely shape id here!
      co_density[sh] = m_imDensityCo.deriv (sh, _C_, sh); // we assume the completely local dependence of \f$\rho\f$ of \f$c\f$
      ConsGravityMethod.template prepare<dim>
        (vConsGravity_c[sh], numSh, vCornerCoords, co_density, m_Gravity);
      co_density[sh] = 0.0;
    }
    for (size_t sh = 0; sh < numSh; sh++)
    { //TODO: use scv, not merely shape id here!
      co_density[sh] = m_imDensityCo.deriv (sh, _T_, sh); // we assume the completely local dependence of \f$\rho\f$ of \f$T\f$
      ConsGravityMethod.template prepare<dim>
        (vConsGravity_T[sh], numSh, vCornerCoords, co_density, m_Gravity);
      co_density[sh] = 0.0;
    }
  }
  UG_CATCH_THROW ("FractTHF_FV1::bulk_ass_JA_elem: Cannot prepare Consistent Gravity or its derivatives.");
 
//  Get the corner pressure
  number vPressure [numCorners];
  for (size_t sh = 0; sh < numSh; sh++)
    vPressure[sh] = u(_P_, sh);
 
//  Compute Diffusion - Dispersion
  MatScale(Diffusion[0], m_imPorosity[0], m_imDiffusion[0]);
  for (size_t ip = 1; ip < numScvf; ip++)
    Diffusion[ip] = Diffusion[0];
  //TODO add Dispersion
 
//  Compute the thermal conductivity
  ThermCond[0] = m_imThermalConductivity[0];
  for (size_t ip = 1; ip < numScvf; ip++)
    ThermCond[ip] = ThermCond[0];
  
//  Compute the Darcy velocity at the ips and the convection shapes
  for (size_t ip = 0; ip < numScvf; ip++)
  {
      this->template compute_ip_Darcy_velocity<TBulkFVGeom, TConsGravity>
          (Vel[ip], ip, bulk_geo, ConsGravityMethod, vConsGravity, vPressure,
            m_imViscosityIP[ip]);
      MatVecMult (DarcyVel[ip], m_imPermeability[0], Vel[ip]);
  }
  
//  compute upwind shapes for transport equation
  if(!m_spUpwind->update(&bulk_geo, DarcyVel, Diffusion, true))
    UG_THROW("FractTHF_FV1::bulk_ass_dA_elem: Cannot compute convection shapes for the transport equation.");
 
//  compute upwind shapes for transport equation
  if(!m_spUpwindEnergy->update(&bulk_geo, DarcyVel, ThermCond, true))
    UG_THROW("FractTHF_FV1::bulk_ass_dA_elem: Cannot compute convection shapes for the energy equation.");
 
//  get a const (!!) reference to the upwind
  const IConvectionShapes<dim>& convShape
    = *const_cast<const IConvectionShapes<dim>*>(m_spUpwind.get());
  
//  get a const (!!) reference to the upwind
  const IConvectionShapes<dim>& convShapeT
    = *const_cast<const IConvectionShapes<dim>*>(m_spUpwindEnergy.get());
  
//  Loop Sub Control Volume Faces (SCVF)
  for (size_t ip = 0; ip < numScvf; ip++)
  {
  //  Get current SCVF
    const typename TBulkFVGeom::SCVF& scvf = bulk_geo.scvf(ip);
 
    // The ip-concentration and the gradient of c at the ip:
    VecSet (grad_c_ip, 0.0); VecSet (grad_T_ip, 0.0);
    for (size_t sh = 0; sh < scvf.num_sh(); sh++)
    {
      VecScaleAppend (grad_c_ip, u(_C_, sh), scvf.global_grad(sh));
      VecScaleAppend (grad_T_ip, u(_T_, sh), scvf.global_grad(sh));
    }
 
  //  Compute the derivative of the viscosity
    if (m_imViscosityIP.constant ())
    {
      pViscosity_c = NULL;
      pViscosity_T = NULL;
    }
    else
    {
      for (size_t sh = 0; sh < numSh; sh++)
      {
        vViscosity_c[sh] = m_imViscosityIP.deriv(ip, _C_, sh);
        vViscosity_T[sh] = m_imViscosityIP.deriv(ip, _T_, sh);
      }
      pViscosity_c = vViscosity_c;
      pViscosity_T = vViscosity_T;
    }
 
  // Compute derivatives of the Darcy velocity
      this->template compute_J_ip_Darcy_velocity <TBulkFVGeom, TConsGravity, numCorners>
          (Vel[ip], Vel_c, Vel_p, Vel_T, ip, bulk_geo, ConsGravityMethod,
            vConsGravity_c, vConsGravity_T, vPressure,
            m_imViscosityIP[ip], pViscosity_c, pViscosity_T);
      for (size_t sh = 0; sh < numSh; sh++)
      {
        MatVecMult (vDDarcyVel_c[sh], m_imPermeability[0], Vel_c[sh]);
        MatVecMult (vDDarcyVel_p[sh], m_imPermeability[0], Vel_p[sh]);
        MatVecMult (vDDarcyVel_T[sh], m_imPermeability[0], Vel_T[sh]);
      }
 
  // Transport Equation
 
   // Loop Shape Functions
    for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
    {
    //  Compute Derivative of Convective Flux
      vDFlux_c[sh] = convShape(ip, sh);
      vDFlux_p[sh] = 0.0;
      vDFlux_T[sh] = 0.0;
 
    //  Derivative w.r.t. Velocity
      for(size_t sh1 = 0; sh1 < scvf.num_sh(); ++sh1)
      {
        vDFlux_c[sh] += u(_C_, sh1) * VecDot(convShape.D_vel(ip, sh1), vDDarcyVel_c[sh]);
        vDFlux_p[sh] += u(_C_, sh1) * VecDot(convShape.D_vel(ip, sh1), vDDarcyVel_p[sh]);
        vDFlux_T[sh] += u(_C_, sh1) * VecDot(convShape.D_vel(ip, sh1), vDDarcyVel_T[sh]);
      }
 
    //TODO: add derivative of Dispersion
 
    //  Add Derivative of Diffusive Flux
      MatVecMult(Dgrad, Diffusion[ip], scvf.global_grad(sh));
      vDFlux_c[sh] -= VecDot(Dgrad, scvf.normal());
    }
 
  //  Handle density in case of full equation
    if(!m_BoussinesqTransport)
    {
    //  Convective Flux
      number flux = 0;
      for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
        flux += convShape(ip, sh) * u(_C_, sh);
 
    //  Diffusive Flux
      MatVecMult(Dgrad, Diffusion[ip], grad_c_ip);
      flux -= VecDot(Dgrad, scvf.normal());
 
    //  Derivative of product
      for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
      {
        vDFlux_c[sh] = m_imDensityIP[ip] * vDFlux_c[sh] +
                 m_imDensityIP.deriv(ip, _C_, sh) * flux;
        vDFlux_p[sh] *= m_imDensityIP[ip];
        vDFlux_T[sh] = m_imDensityIP[ip] * vDFlux_T[sh] +
                 m_imDensityIP.deriv(ip, _T_, sh) * flux;
      }
    }
 
  //  Add Flux contribution
    for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
    {
      J(_C_, scvf.from(), _C_, sh) += vDFlux_c[sh];
      J(_C_, scvf.to(),   _C_, sh) -= vDFlux_c[sh];
      J(_C_, scvf.from(), _P_, sh) += vDFlux_p[sh];
      J(_C_, scvf.to(),   _P_, sh) -= vDFlux_p[sh];
      J(_C_, scvf.from(), _T_, sh) += vDFlux_T[sh];
      J(_C_, scvf.to(),   _T_, sh) -= vDFlux_T[sh];
    }
 
 
  // Flow Equation
 
    for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
    {
      vDFlux_c[sh] = VecDot(vDDarcyVel_c[sh], scvf.normal());
      vDFlux_p[sh] = VecDot(vDDarcyVel_p[sh], scvf.normal());
      vDFlux_T[sh] = VecDot(vDDarcyVel_T[sh], scvf.normal());
    }
 
    number flux = VecDot(DarcyVel[ip], scvf.normal());
    
    if(!m_BoussinesqFlow)
    {
      for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
      {
        vDFlux_c[sh] = m_imDensityIP[ip] * vDFlux_c[sh] +
                 m_imDensityIP.deriv(ip, _C_, sh) * flux;
        vDFlux_p[sh] *= m_imDensityIP[ip];
        vDFlux_T[sh] = m_imDensityIP[ip] * vDFlux_T[sh] +
                 m_imDensityIP.deriv(ip, _T_, sh) * flux;
      }
      
      flux *= m_imDensityIP[ip];
    }
 
  //  Add Flux contribution
    for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
    {
      J(_P_, scvf.from(), _C_, sh) += vDFlux_c[sh];
      J(_P_, scvf.to(),   _C_, sh) -= vDFlux_c[sh];
      J(_P_, scvf.from(), _P_, sh) += vDFlux_p[sh];
      J(_P_, scvf.to(),   _P_, sh) -= vDFlux_p[sh];
      J(_P_, scvf.from(), _T_, sh) += vDFlux_T[sh];
      J(_P_, scvf.to(),   _T_, sh) -= vDFlux_T[sh];
    }
  
  //  Volume stabilization: Add the continuity equation to the transport equation
  //  REMARK: We do not compute derivatives of the defect of the conf. eq. here because this defect tends to 0.
    if(m_bVolStabDataActive)
    {
      m_spVolStabData->stiff(pElem->vertex(scvf.from())) -= flux;
      m_spVolStabData->stiff(pElem->vertex(scvf.to())) += flux;
    }
 
 
  // Energy Equation
 
   // Loop Shape Functions
    for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
    {
    //  Compute Derivative of Convective Flux
      vDFlux_c[sh] = 0.0;
      vDFlux_p[sh] = 0.0;
      vDFlux_T[sh] = convShape(ip, sh);
 
    //  Derivative w.r.t. Velocity
      for(size_t sh1 = 0; sh1 < scvf.num_sh(); ++sh1)
      {
        vDFlux_c[sh] += u(_T_, sh1) * VecDot(convShapeT.D_vel(ip, sh1), vDDarcyVel_c[sh]);
        vDFlux_p[sh] += u(_T_, sh1) * VecDot(convShapeT.D_vel(ip, sh1), vDDarcyVel_p[sh]);
        vDFlux_T[sh] += u(_T_, sh1) * VecDot(convShapeT.D_vel(ip, sh1), vDDarcyVel_T[sh]);
      }
 
      vDFlux_c[sh] *= m_imHeatCapacityFluid;
      vDFlux_p[sh] *= m_imHeatCapacityFluid;
      vDFlux_T[sh] *= m_imHeatCapacityFluid;
    }
 
  //  Handle density in case of full equation
    if(!m_BoussinesqEnergy)
    {
    //  Convective Flux
      number flux = 0;
      for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
        flux += convShapeT(ip, sh) * u(_T_, sh);
      flux *= m_imHeatCapacityFluid;
 
    //  Derivative of product
      for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
      {
        vDFlux_c[sh] = m_imDensityIP[ip] * vDFlux_c[sh] +
                 m_imDensityIP.deriv(ip, _C_, sh) * flux;
        vDFlux_p[sh] *= m_imDensityIP[ip];
        vDFlux_T[sh] = m_imDensityIP[ip] * vDFlux_T[sh] +
                 m_imDensityIP.deriv(ip, _T_, sh) * flux;
      }
    }
    else
    {
      for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
      {
        vDFlux_c[sh] *= m_BoussinesqDensity;
        vDFlux_p[sh] *= m_BoussinesqDensity;
        vDFlux_T[sh] *= m_BoussinesqDensity;
      }
    }
 
  //  Loop Shape Functions
    for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
    {
    //  Add Derivative of Diffusive Flux
      MatVecMult(Dgrad, ThermCond[ip], scvf.global_grad(sh));
      vDFlux_T[sh] -= VecDot(Dgrad, scvf.normal());
    }
 
  //  Add Flux contribution
    for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
    {
      J(_T_, scvf.from(), _C_, sh) += vDFlux_c[sh];
      J(_T_, scvf.to(),   _C_, sh) -= vDFlux_c[sh];
      J(_T_, scvf.from(), _P_, sh) += vDFlux_p[sh];
      J(_T_, scvf.to(),   _P_, sh) -= vDFlux_p[sh];
      J(_T_, scvf.from(), _T_, sh) += vDFlux_T[sh];
      J(_T_, scvf.to(),   _T_, sh) -= vDFlux_T[sh];
    }
  }
  
//  Assemble the singular sources and sinks
    if (m_sss_mngr.valid () && (m_sss_mngr->num_points () != 0 || m_sss_mngr->num_lines () != 0))
    {
      typedef typename domain_type::position_accessor_type t_pos_accessor;
      typedef typename THSingularSourcesAndSinks<dim>::template
        point_iterator<TElem,t_pos_accessor,TBulkFVGeom> t_pnt_sss_iter;
      typedef typename THSingularSourcesAndSinks<dim>::template
        line_iterator<TElem,t_pos_accessor,TBulkFVGeom> t_lin_sss_iter;
      
    t_pos_accessor& aaPos = this->domain()->position_accessor ();
    Grid& grid = (Grid&) *this->domain()->grid ();
    
    for(size_t ip = 0; ip < bulk_geo.num_scv(); ip++)
    {
      size_t co = bulk_geo.scv(ip).node_id ();
      
    //  point sources
      for (t_pnt_sss_iter pnt (m_sss_mngr.get (), pElem, grid, aaPos, bulk_geo, co);
        ! pnt.is_over (); ++pnt)
      {
        FVPointSourceOrSink<dim, th_point_sss_data<dim> > * pnt_sss = *pnt;
        if (! pnt_sss->marked_for (pElem, co))
          continue;
        pnt_sss->compute (pnt_sss->position (), this->time (), -1); //TODO: set the subset id instead of -1
        ass_sss_JA_elem (J, u, pElem, co, pnt_sss->intensity (), pnt_sss->concentration (), pnt_sss->temperature ());
      }
      
    //  line sources
      for (t_lin_sss_iter line (m_sss_mngr.get (), pElem, grid, aaPos, bulk_geo, co);
        ! line.is_over (); ++line)
      {
        FVLineSourceOrSink<dim, th_line_sss_data<dim> > * line_sss = *line;
        number len = VecDistance (line.seg_start (), line.seg_end ());
        line_sss->compute (line.seg_start (), this->time (), -1); //TODO: set the subset id instead of -1
        ass_sss_JA_elem (J, u, pElem, co, line_sss->intensity (), line_sss->concentration (), line_sss->temperature (), len);
      }
    }
    }
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::fract_ass_JA_elem
(
  LocalMatrix & J, 
  const LocalVector & u, 
  TElem * pElem, 
  const position_type vCornerCoords [] 
)
{
  number half_fr_width = m_imAperture[0] / 2;
  
  const size_t numSh = m_pFractGeo->num_sh();
  const size_t numScvf = m_pFractGeo->num_scvf();
 
  MathMatrix<dim,dim> Diffusion[TFractFVGeom::maxNumSCVF];
  MathMatrix<dim,dim> ThermCond[TFractFVGeom::maxNumSCVF];
  MathVector<dim> Dgrad;
  MathVector<dim> grad_c_ip, grad_T_ip;
  number vDFlux_c [maxFractSideCorners], vDFlux_p [maxFractSideCorners], vDFlux_T [maxFractSideCorners];
 
//  Get the corner values
  number vConcentration [maxFractSideCorners];
  number vTemperature [maxFractSideCorners];
  for (size_t sh = 0; sh < numSh; sh++)
  {
    size_t co = m_innerSideCo[sh]; // co: index of the corner in the (degenerated) full-dim. element
    vConcentration[sh] = u(_C_, co);
    vTemperature[sh] = u(_T_, co);
  }
  
//  Compute Diffusion - Dispersion
  MatSet (Diffusion[0], 0);
  MatDiagSet (Diffusion[0], m_imPorosity[0] * m_imFractDiffusion[0]);
  for (size_t ip = 1; ip < numScvf; ip++)
    Diffusion[ip] = Diffusion[0];
  //TODO add Dispersion
  
//  Compute the thermal conductivity
  MatSet (ThermCond[0], 0);
  MatDiagSet (ThermCond[0], m_imFractThermCond[0]);
  for (size_t ip = 1; ip < numScvf; ip++)
    ThermCond[ip] = ThermCond[0];
  
//  compute upwind shapes for transport equation
  if(!m_spUpwind->update(m_pFractGeo, m_imFractDarcyVelIP.values(), Diffusion, true))
    UG_THROW("FractTHF_FV1::fract_ass_dA_elem: Cannot compute convection shapes for the transport equation.");
 
//  compute upwind shapes for transport equation
  if(!m_spUpwindEnergy->update(m_pFractGeo, m_imFractDarcyVelIP.values(), ThermCond, true))
    UG_THROW("FractTHF_FV1::fract_ass_dA_elem: Cannot compute convection shapes for the energy equation.");
 
//  get a const (!!) reference to the upwind
  const IConvectionShapes<dim>& convShape
    = *const_cast<const IConvectionShapes<dim>*>(m_spUpwind.get());
  
//  get a const (!!) reference to the upwind
  const IConvectionShapes<dim>& convShapeT
    = *const_cast<const IConvectionShapes<dim>*>(m_spUpwindEnergy.get());
  
//  Loop Sub Control Volume Faces (SCVF)
  for (size_t ip = 0; ip < numScvf; ip++)
  {
  //  Get current SCVF
    const typename TFractFVGeom::SCVF& scvf = m_pFractGeo->scvf(ip);
 
    // The ip-concentration and the gradient of c at the ip:
    VecSet (grad_c_ip, 0.0); VecSet (grad_T_ip, 0.0);
    for (size_t sh = 0; sh < scvf.num_sh(); sh++)
    {
      VecScaleAppend (grad_c_ip, vConcentration[sh], scvf.global_grad(sh));
      VecScaleAppend (grad_T_ip, vTemperature[sh], scvf.global_grad(sh));
    }
 
    const MathVector<dim>* vDDarcyVel_c = m_imFractDarcyVelIP.deriv(ip, _C_);
    const MathVector<dim>* vDDarcyVel_p = m_imFractDarcyVelIP.deriv(ip, _P_);
    const MathVector<dim>* vDDarcyVel_T = m_imFractDarcyVelIP.deriv(ip, _T_);
 
  // Transport Equation
 
   // Loop Shape Functions
    for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
    {
      size_t co = m_innerSideCo[sh]; // co: index of the corner in the (degenerated) full-dim. element
      
    //  Compute Derivative of Convective Flux
      vDFlux_c[sh] = convShape(ip, sh);
      vDFlux_p[sh] = 0.0;
      vDFlux_T[sh] = 0.0;
 
    //  Derivative w.r.t. Velocity
      for(size_t sh1 = 0; sh1 < scvf.num_sh(); ++sh1)
      {
        vDFlux_c[sh] += vConcentration[sh1] * VecDot(convShape.D_vel(ip, sh1), vDDarcyVel_c[co]);
        vDFlux_p[sh] += vConcentration[sh1] * VecDot(convShape.D_vel(ip, sh1), vDDarcyVel_p[co]);
        vDFlux_T[sh] += vConcentration[sh1] * VecDot(convShape.D_vel(ip, sh1), vDDarcyVel_T[co]);
      }
    //TODO: add derivative of Dispersion
 
    //  Add Derivative of Diffusive Flux
      MatVecMult(Dgrad, Diffusion[ip], scvf.global_grad(sh));
      vDFlux_c[sh] -= VecDot(Dgrad, scvf.normal());
    }
    
  //  Handle density in case of full equation
    if(!m_BoussinesqTransport)
    {
    //  Convective Flux
      number flux = 0;
      for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
        flux += convShape(ip, sh) * vConcentration[sh];
 
    //  Diffusive Flux
      MatVecMult(Dgrad, Diffusion[ip], grad_c_ip);
      flux -= VecDot(Dgrad, scvf.normal());
 
    //  Derivative of product
      for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
      {
        vDFlux_c[sh] = m_imDensityIP[ip] * vDFlux_c[sh] +
                (m_imDensityIP.constant() ? 0.0 : m_imDensityIP.deriv(ip, _C_, m_innerSideCo[sh])) * flux;
        vDFlux_p[sh] *= m_imDensityIP[ip];
        vDFlux_T[sh] = m_imDensityIP[ip] * vDFlux_T[sh] +
                (m_imDensityIP.constant() ? 0.0 : m_imDensityIP.deriv(ip, _T_, m_innerSideCo[sh])) * flux;
      }
    }
 
  //  Multiply by the fracture width
    for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
    {
      vDFlux_c[sh] *= half_fr_width; vDFlux_p[sh] *= half_fr_width; vDFlux_T[sh] *= half_fr_width;
    }
    
  //  Add Flux contribution
    for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
    {
      size_t co = m_innerSideCo[sh];
      size_t co_from = m_innerSideCo[scvf.from()];
      size_t co_to = m_innerSideCo[scvf.to()];
      
      J(_C_, co_from, _C_, co) += vDFlux_c[sh];
      J(_C_, co_to,   _C_, co) -= vDFlux_c[sh];
      J(_C_, co_from, _P_, co) += vDFlux_p[sh];
      J(_C_, co_to,   _P_, co) -= vDFlux_p[sh];
      J(_C_, co_from, _T_, co) += vDFlux_T[sh];
      J(_C_, co_to,   _T_, co) -= vDFlux_T[sh];
    }
 
  // Flow Equation
 
    for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
    {
      size_t co = m_innerSideCo[sh];
      vDFlux_c[sh] = VecDot(vDDarcyVel_c[co], scvf.normal());
      vDFlux_p[sh] = VecDot(vDDarcyVel_p[co], scvf.normal());
      vDFlux_T[sh] = VecDot(vDDarcyVel_T[co], scvf.normal());
    }
 
    number flux = VecDot(m_imFractDarcyVelIP[ip], scvf.normal());
    
    if(!m_BoussinesqFlow)
    {
      for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
      {
        vDFlux_c[sh] = m_imDensityIP[ip] * vDFlux_c[sh] +
                (m_imDensityIP.constant() ? 0.0 : m_imDensityIP.deriv(ip, _C_, m_innerSideCo[sh])) * flux;
        vDFlux_p[sh] *= m_imDensityIP[ip];
        vDFlux_T[sh] = m_imDensityIP[ip] * vDFlux_T[sh] +
                (m_imDensityIP.constant() ? 0.0 : m_imDensityIP.deriv(ip, _T_, m_innerSideCo[sh])) * flux;
      }
      
      flux *= m_imDensityIP[ip];
    }
 
  //  Multiply by the fracture width
    for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
    {
      vDFlux_c[sh] *= half_fr_width; vDFlux_p[sh] *= half_fr_width; vDFlux_T[sh] *= half_fr_width;
    }
    
  //  Add Flux contribution
    for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
    {
      size_t co = m_innerSideCo[sh];
      size_t co_from = m_innerSideCo[scvf.from()];
      size_t co_to = m_innerSideCo[scvf.to()];
      
      J(_P_, co_from, _C_, co) += vDFlux_c[sh];
      J(_P_, co_to,   _C_, co) -= vDFlux_c[sh];
      J(_P_, co_from, _P_, co) += vDFlux_p[sh];
      J(_P_, co_to,   _P_, co) -= vDFlux_p[sh];
      J(_P_, co_from, _T_, co) += vDFlux_T[sh];
      J(_P_, co_to,   _T_, co) -= vDFlux_T[sh];
    }
  
  //  Volume stabilization: Add the continuity equation to the transport equation
  //  REMARK: We do not compute derivatives of the defect of the conf. eq. here because this defect tends to 0.
    if(m_bVolStabDataActive)
    {
      flux *= half_fr_width;
      m_spVolStabData->stiff(pElem->vertex(m_innerSideCo[scvf.from()])) -= flux;
      m_spVolStabData->stiff(pElem->vertex(m_innerSideCo[scvf.to()])) += flux;
    }
 
 
  // Energy Equation
 
   // Loop Shape Functions
    for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
    {
      size_t co = m_innerSideCo[sh]; // co: index of the corner in the (degenerated) full-dim. element
      
    //  Compute Derivative of Convective Flux
      vDFlux_c[sh] = 0.0;
      vDFlux_p[sh] = 0.0;
      vDFlux_T[sh] = convShape(ip, sh);
 
    //  Derivative w.r.t. Velocity
      for(size_t sh1 = 0; sh1 < scvf.num_sh(); ++sh1)
      {
        vDFlux_c[sh] += vTemperature[sh1] * VecDot(convShapeT.D_vel(ip, sh1), vDDarcyVel_c[co]);
        vDFlux_p[sh] += vTemperature[sh1] * VecDot(convShapeT.D_vel(ip, sh1), vDDarcyVel_p[co]);
        vDFlux_T[sh] += vTemperature[sh1] * VecDot(convShapeT.D_vel(ip, sh1), vDDarcyVel_T[co]);
      }
 
      vDFlux_c[sh] *= m_imHeatCapacityFluid;
      vDFlux_p[sh] *= m_imHeatCapacityFluid;
      vDFlux_T[sh] *= m_imHeatCapacityFluid;
    }
 
  //  Handle density in case of full equation
    if(!m_BoussinesqEnergy)
    {
    //  Convective Flux
      number flux = 0;
      for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
        flux += convShapeT(ip, sh) * vTemperature[sh];
      flux *= m_imHeatCapacityFluid;
 
    //  Derivative of product
      for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
      {
        vDFlux_c[sh] = m_imDensityIP[ip] * vDFlux_c[sh] +
                (m_imDensityIP.constant() ? 0.0 : m_imDensityIP.deriv(ip, _C_, m_innerSideCo[sh])) * flux;
        vDFlux_p[sh] *= m_imDensityIP[ip];
        vDFlux_T[sh] = m_imDensityIP[ip] * vDFlux_T[sh] +
                (m_imDensityIP.constant() ? 0.0 : m_imDensityIP.deriv(ip, _T_, m_innerSideCo[sh])) * flux;
      }
    }
    else
    {
      for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
      {
        vDFlux_c[sh] *= m_BoussinesqDensity;
        vDFlux_p[sh] *= m_BoussinesqDensity;
        vDFlux_T[sh] *= m_BoussinesqDensity;
      }
    }
 
  //  Loop Shape Functions
    for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
    {
    //  Add Derivative of Diffusive Flux
      MatVecMult(Dgrad, ThermCond[ip], scvf.global_grad(sh));
      vDFlux_T[sh] -= VecDot(Dgrad, scvf.normal());
    }
 
  //  Multiply by the fracture width
    for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
    {
      vDFlux_c[sh] *= half_fr_width; vDFlux_p[sh] *= half_fr_width; vDFlux_T[sh] *= half_fr_width;
    }
    
  //  Add Flux contribution
    for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
    {
      size_t co = m_innerSideCo[sh];
      size_t co_from = m_innerSideCo[scvf.from()];
      size_t co_to = m_innerSideCo[scvf.to()];
      
      J(_T_, co_from, _C_, co) += vDFlux_c[sh];
      J(_T_, co_to,   _C_, co) -= vDFlux_c[sh];
      J(_T_, co_from, _P_, co) += vDFlux_p[sh];
      J(_T_, co_to,   _P_, co) -= vDFlux_p[sh];
      J(_T_, co_from, _T_, co) += vDFlux_T[sh];
      J(_T_, co_to,   _T_, co) -= vDFlux_T[sh];
    }
  }
  
//  Assemble the singular sources and sinks
    if (m_sss_mngr.valid () && m_sss_mngr->num_lines () != 0)
    {
      typedef typename domain_type::position_accessor_type t_pos_accessor;
      typedef typename THSingularSourcesAndSinks<dim>::template
        line_iterator<side_type,t_pos_accessor,TFractFVGeom> t_lin_sss_iter;
      
    t_pos_accessor& aaPos = this->domain()->position_accessor ();
    Grid& grid = (Grid&) *this->domain()->grid ();
    
    for(size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
    {
    //  Get the corner of the face
      size_t side_co = m_pFractGeo->scv(ip).node_id ();
    //  Get associated node of the element (not side!)
      size_t co = m_innerSideCo [m_pFractGeo->scv(ip).node_id ()];
      
    //  line sources (that correspond to the point sources)
      for (t_lin_sss_iter line (m_sss_mngr.get (), m_innerFractSide, grid, aaPos, *m_pFractGeo, side_co);
        ! line.is_over (); ++line)
      {
        FVLineSourceOrSink<dim, th_line_sss_data<dim> > * line_sss = *line;
        if (! line_sss->marked_for (m_innerFractSide, side_co))
          continue;
        line_sss->compute (line.seg_start (), this->time (), -1); //TODO: set the subset id instead of -1
        ass_sss_JA_elem (J, u, pElem, co, line_sss->intensity () / 2, line_sss->concentration (), line_sss->temperature (), 0.5);
        /* Remark: "0.5" because the source is taken into account twise. */
      }
    }
  }
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::fract_bulk_ass_JA_elem
(
  LocalMatrix & J, 
  const LocalVector & u, 
  TElem * pElem, 
  const position_type vCornerCoords [] 
)
{
  const number orthDiffusion = m_imPorosity[0] * m_imOrthoDiffusion[0];
  
  const number half_fr_width = m_imAperture[0] / 2;
  
//  loop over the corners of the inner side
  for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
  {
  //  Get current SCV
    const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv(ip);
    number s = scv.volume ();
 
  //  Get associated node of the element (not side!)
    const int co = m_innerSideCo [scv.node_id()];
    
  //  Get the corner values
    const number orthC_f = u(_C_, co);
    const number orthC_m = u(_C_, m_assCo[co]);
    const number orthP_f = u(_P_, co);
    const number orthP_m = u(_P_, m_assCo[co]);
    const number orthT_f = u(_T_, co);
    const number orthT_m = u(_T_, m_assCo[co]);
    
    const number fractDensity = m_imDensityCo[co];
    const number orthDensity = m_imDensityCo [m_assCo[co]];
    number D_fractDensity_c, D_orthDensity_c, D_fractDensity_T, D_orthDensity_T;
    if (! m_imDensityCo.constant())
    {
      D_fractDensity_c = m_imDensityCo.deriv (co, _C_, co);
      D_fractDensity_T = m_imDensityCo.deriv (co, _T_, co);
      D_orthDensity_c = m_imDensityCo.deriv (m_assCo[co], _C_, m_assCo[co]);
      D_orthDensity_T = m_imDensityCo.deriv (m_assCo[co], _T_, m_assCo[co]);
    }
    else
      D_fractDensity_c = D_fractDensity_T = D_orthDensity_c = D_orthDensity_T = 0;
  
    const number orthViscosity = m_imViscosityCo [m_assCo[co]];
    number D_orthViscosity_c, D_orthViscosity_T;
    if (! m_imViscosityCo.constant())
    {
      D_orthViscosity_c = m_imViscosityCo.deriv (m_assCo[co], _C_, m_assCo[co]);
      D_orthViscosity_T = m_imViscosityCo.deriv (m_assCo[co], _T_, m_assCo[co]);
    }
    else
      D_orthViscosity_c = D_orthViscosity_T = 0;
    
  //  Flux velocity between the fracture and the bulk medium
    
    number orthVelocity = (m_orthGravity * (orthDensity - fractDensity)
                    - (orthP_m - orthP_f) / half_fr_width)
                      * m_imOrthoPermeability[0];
 
    number D_orthVelocity [3], D_orthVelocity_fr [3]; // 3 components: w.r.t. _C_, _P_ and _T_ resp.
    
        D_orthVelocity [_C_] = (m_orthGravity * D_orthDensity_c
          * m_imOrthoPermeability[0] * orthViscosity
            - orthVelocity * D_orthViscosity_c) / orthViscosity / orthViscosity;
        
        D_orthVelocity_fr [_C_] = - m_orthGravity * D_fractDensity_c
          * m_imOrthoPermeability[0] / orthViscosity;
        
        D_orthVelocity [_P_] = - m_imOrthoPermeability[0]
          / half_fr_width / orthViscosity;
        
        D_orthVelocity_fr [_P_] = - D_orthVelocity [_P_];
        
        D_orthVelocity [_T_] = (m_orthGravity * D_orthDensity_T
          * m_imOrthoPermeability[0] * orthViscosity
            - orthVelocity * D_orthViscosity_T) / orthViscosity / orthViscosity;
        
        D_orthVelocity_fr [_T_] = - m_orthGravity * D_fractDensity_T
          * m_imOrthoPermeability[0] / orthViscosity;
        
        orthVelocity /= orthViscosity;
    
      number flux, D_flux [3], D_flux_fr [3]; // 3 components of the derivatives: w.r.t. _C_, _P_ and _T_ resp.
 
  // Transport Equation
 
    D_flux [_C_] = orthDiffusion / half_fr_width;
    D_flux_fr [_C_] = - D_flux [_C_];
    D_flux [_T_] = D_flux_fr [_T_] = 0;
    /* We use the full upwind here: */
    if (orthVelocity >= 0)
    {
      D_flux [_C_] -= D_orthVelocity [_C_] * orthC_f;
      D_flux_fr [_C_] -= orthVelocity + D_orthVelocity_fr [_C_] * orthC_f;
      D_flux [_P_] = - D_orthVelocity [_P_] * orthC_f;
      D_flux_fr [_P_] = - D_orthVelocity_fr [_P_] * orthC_f;
      D_flux [_T_] = - D_orthVelocity [_T_] * orthC_f;
      D_flux_fr [_T_] = - D_orthVelocity_fr [_T_] * orthC_f;
    }
    else
    {
      D_flux [_C_] -= orthVelocity + D_orthVelocity [_C_] * orthC_m;
      D_flux_fr [_C_] -= D_orthVelocity_fr [_C_] * orthC_m;
      D_flux [_P_] = - D_orthVelocity [_P_] * orthC_m;
      D_flux_fr [_P_] = - D_orthVelocity_fr [_P_] * orthC_m;
      D_flux [_T_] = - D_orthVelocity [_T_] * orthC_m;
      D_flux_fr [_T_] = - D_orthVelocity_fr [_T_] * orthC_m;
    }
    
    if (! m_BoussinesqTransport)
    {
      flux = orthDiffusion * (orthC_m - orthC_f) / half_fr_width;
      /* We use the full upwind here: */
      flux -= orthVelocity * ((orthVelocity >= 0)? orthC_f : orthC_m);
 
      D_flux [_C_] = D_flux [_C_] * orthDensity + flux * D_orthDensity_c;
      D_flux_fr [_C_] *= orthDensity;
      D_flux [_P_] *= orthDensity;
      D_flux_fr [_P_] *= orthDensity;
      D_flux [_T_] = D_flux [_T_] * orthDensity + flux * D_orthDensity_T;
      D_flux_fr [_T_] *= orthDensity;
    }
    
    J(_C_, m_assCo [co], _C_, m_assCo [co]) += D_flux [_C_] * s;
    J(_C_, m_assCo [co], _P_, m_assCo [co]) += D_flux [_P_] * s;
    J(_C_, m_assCo [co], _T_, m_assCo [co]) += D_flux [_T_] * s;
    J(_C_, m_assCo [co], _C_, co) += D_flux_fr [_C_] * s;
    J(_C_, m_assCo [co], _P_, co) += D_flux_fr [_P_] * s;
    J(_C_, m_assCo [co], _T_, co) += D_flux_fr [_T_] * s;
    
    J(_C_, co, _C_, m_assCo [co]) -= D_flux [_C_] * s;
    J(_C_, co, _P_, m_assCo [co]) -= D_flux [_P_] * s;
    J(_C_, co, _T_, m_assCo [co]) -= D_flux [_T_] * s;
    J(_C_, co, _C_, co) -= D_flux_fr [_C_] * s;
    J(_C_, co, _P_, co) -= D_flux_fr [_P_] * s;
    J(_C_, co, _T_, co) -= D_flux_fr [_T_] * s;
 
  // Flow Equation
 
    if (! m_BoussinesqFlow)
    {
      D_flux [_C_] = (D_orthVelocity [_C_] * orthDensity
       + orthVelocity * D_orthDensity_c) * s;
      D_flux_fr [_C_] = D_orthVelocity_fr [_C_] * orthDensity * s;
    }
    else
    {
      D_flux [_C_] = D_orthVelocity [_C_] * s;
      D_flux_fr [_C_] = D_orthVelocity_fr [_C_] * s;
    }
    J(_P_, m_assCo [co], _C_, m_assCo [co]) -= D_flux [_C_];
    J(_P_, co, _C_, m_assCo [co]) += D_flux [_C_];
    J(_P_, m_assCo [co], _C_, co) -= D_flux_fr [_C_];
    J(_P_, co, _C_, co) += D_flux_fr [_C_];
    
    number rho_x_s = s;
    if (! m_BoussinesqFlow) rho_x_s *= orthDensity;
    
    flux = D_orthVelocity [_P_] * rho_x_s;
    J(_P_, m_assCo [co], _P_, m_assCo [co]) -= flux;
    J(_P_, co, _P_, m_assCo [co]) += flux;
    
    flux = D_orthVelocity_fr [_P_] * rho_x_s;
    J(_P_, m_assCo [co], _P_, co) -= flux;
    J(_P_, co, _P_, co) += flux;
 
    if (! m_BoussinesqFlow)
    {
      D_flux [_T_] = (D_orthVelocity [_T_] * orthDensity
       + orthVelocity * D_orthDensity_T) * s;
      D_flux_fr [_T_] = D_orthVelocity_fr [_T_] * orthDensity * s;
    }
    else
    {
      D_flux [_T_] = D_orthVelocity [_T_] * s;
      D_flux_fr [_T_] = D_orthVelocity_fr [_T_] * s;
    }
    J(_P_, m_assCo [co], _T_, m_assCo [co]) -= D_flux [_T_];
    J(_P_, co, _T_, m_assCo [co]) += D_flux [_T_];
    J(_P_, m_assCo [co], _T_, co) -= D_flux_fr [_T_];
    J(_P_, co, _T_, co) += D_flux_fr [_T_];
  
  //  Volume stabilization: Add the continuity equation to the transport equation
  //  REMARK: We do not compute derivatives of the defect of the conf. eq. here because this defect tends to 0.
    if(m_bVolStabDataActive)
    {
      flux = orthVelocity * rho_x_s;
      m_spVolStabData->stiff(pElem->vertex(m_assCo[co])) += flux;
      m_spVolStabData->stiff(pElem->vertex(co)) -= flux;
    }
 
  // Energy Equation
 
    /* We use the full upwind here: */
    if (orthVelocity >= 0)
    {
      D_flux [_C_] = - D_orthVelocity [_C_] * orthT_f * m_imHeatCapacityFluid;
      D_flux_fr [_C_] = - D_orthVelocity_fr [_C_] * orthT_f * m_imHeatCapacityFluid;
      D_flux [_P_] = - D_orthVelocity [_P_] * orthT_f * m_imHeatCapacityFluid;
      D_flux_fr [_P_] = - D_orthVelocity_fr [_P_] * orthT_f * m_imHeatCapacityFluid;
      D_flux [_T_] = D_orthVelocity [_T_] * orthT_f * m_imHeatCapacityFluid;
      D_flux_fr [_T_] = orthVelocity + D_orthVelocity_fr [_T_] * orthT_f * m_imHeatCapacityFluid;
    }
    else
    {
      D_flux [_C_] = - D_orthVelocity [_C_] * orthT_m * m_imHeatCapacityFluid;
      D_flux_fr [_C_] = - D_orthVelocity_fr [_C_] * orthT_m * m_imHeatCapacityFluid;
      D_flux [_P_] = - D_orthVelocity [_P_] * orthT_m * m_imHeatCapacityFluid;
      D_flux_fr [_P_] = - D_orthVelocity_fr [_P_] * orthT_m * m_imHeatCapacityFluid;
      D_flux [_T_] = orthVelocity + D_orthVelocity [_T_] * orthT_m * m_imHeatCapacityFluid;
      D_flux_fr [_T_] = D_orthVelocity_fr [_T_] * orthT_m * m_imHeatCapacityFluid;
    }
    
    if (! m_BoussinesqEnergy)
    {
      /* We use the full upwind here: */
      flux = - orthVelocity * ((orthVelocity >= 0)? orthT_f : orthT_m)
          * m_imHeatCapacityFluid * orthDensity;
      
      D_flux [_C_] = D_flux [_C_] * orthDensity + flux * D_orthDensity_c;
      D_flux_fr [_C_] *= orthDensity;
      D_flux [_P_] *= orthDensity;
      D_flux_fr [_P_] *= orthDensity;
      D_flux [_T_] = D_flux [_T_] * orthDensity + flux * D_orthDensity_T;
      D_flux_fr [_T_] *= orthDensity;
    }
    else
    {
      D_flux [_C_] *= m_BoussinesqDensity;
      D_flux_fr [_C_] *= m_BoussinesqDensity;
      D_flux [_P_] *= m_BoussinesqDensity;
      D_flux_fr [_P_] *= m_BoussinesqDensity;
      D_flux [_T_] *= m_BoussinesqDensity;
      D_flux_fr [_T_] *= m_BoussinesqDensity;
    }
    
    D_flux [_T_] += m_imOrthoThermCond[0] / half_fr_width;
    D_flux_fr [_T_] -= m_imOrthoThermCond[0] / half_fr_width;
    
    J(_T_, m_assCo [co], _C_, m_assCo [co]) += D_flux [_C_] * s;
    J(_T_, m_assCo [co], _P_, m_assCo [co]) += D_flux [_P_] * s;
    J(_T_, m_assCo [co], _T_, m_assCo [co]) += D_flux [_T_] * s;
    J(_T_, m_assCo [co], _C_, co) += D_flux_fr [_C_] * s;
    J(_T_, m_assCo [co], _P_, co) += D_flux_fr [_P_] * s;
    J(_T_, m_assCo [co], _T_, co) += D_flux_fr [_T_] * s;
    
    J(_T_, co, _C_, m_assCo [co]) -= D_flux [_C_] * s;
    J(_T_, co, _P_, m_assCo [co]) -= D_flux [_P_] * s;
    J(_T_, co, _T_, m_assCo [co]) -= D_flux [_T_] * s;
    J(_T_, co, _C_, co) -= D_flux_fr [_C_] * s;
    J(_T_, co, _P_, co) -= D_flux_fr [_P_] * s;
    J(_T_, co, _T_, co) -= D_flux_fr [_T_] * s;
  }
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::ass_sss_dA_elem
(
  LocalVector& d, 
  const LocalVector& u, 
  TElem * pElem, 
  size_t co, 
  number intensity, 
  number concentration, 
  number temperature, 
  number factor 
)
{
  if (intensity != 0)
  { // fluid source/sink
    intensity *= factor;
    if (intensity > 0)
    {
    // fluid source
      number density = 1000 + 000 * concentration; //FIXME!!
      d(_C_, co) -= m_BoussinesqTransport? intensity * concentration : intensity * density * concentration;
      d(_P_, co) -= m_BoussinesqFlow? intensity : intensity * density;
      d(_T_, co) -= intensity * (m_BoussinesqEnergy? m_BoussinesqDensity : density) * m_imHeatCapacityFluid * temperature;
      if(m_bVolStabDataActive)
        m_spVolStabData->stiff(pElem->vertex(co)) -= m_BoussinesqFlow? intensity : intensity * density;
    }
    else
    {
    // fluid sink
      d(_C_, co) -= m_BoussinesqTransport? intensity * u(_C_, co) : intensity * m_imDensityCo[co] * u(_C_, co);
      d(_P_, co) -= m_BoussinesqFlow? intensity : intensity * m_imDensityCo[co];
      d(_T_, co) -= intensity * (m_BoussinesqEnergy? m_BoussinesqDensity : m_imDensityCo[co]) * m_imHeatCapacityFluid * u(_T_, co);
      if(m_bVolStabDataActive)
        m_spVolStabData->stiff(pElem->vertex(co))
        -= m_BoussinesqFlow? intensity : intensity * m_imDensityCo[co];
    }
  }
  else
  { // heat source
    if (temperature < 0)
      UG_THROW ("FractTHF_FV1: Heat sinks are not implemented.");
    d(_T_, co) -= temperature * factor;
  }
}
  
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::ass_sss_JA_elem
(
  LocalMatrix& J, 
  const LocalVector& u, 
  TElem * pElem, 
  size_t co, 
  number intensity, 
  number concentration, 
  number temperature, 
  number factor 
)
{
  if (intensity != 0)
  {
    intensity *= factor;
    if (intensity > 0)
    {
    // fluid source: no contribution to the maxrix
      number density  = 1000 + 000 * concentration; //FIXME!!
      if(m_bVolStabDataActive)
        m_spVolStabData->stiff(pElem->vertex(co))
          -= m_BoussinesqFlow? intensity : intensity * density;
    }
    else
    {
    // fluid sink: add the intensity to the matrix
      if (! m_BoussinesqTransport)
      {
        if (!m_imDensityCo.constant ())
        {
          J(_C_, co, _C_, co)
            -= intensity * (m_imDensityCo[co] + m_imDensityCo.deriv (co, _C_, co) * u (_C_, co));
          J(_C_, co, _T_, co) -= intensity * m_imDensityCo.deriv(co, _T_, co) * m_imDensityCo[co];
        }
        else
          J(_C_, co, _C_, co) -= intensity * m_imDensityCo[co];
      }
      else
        J(_C_, co, _C_, co) -= intensity;
    
      if (! m_BoussinesqFlow && ! m_imDensityCo.constant ())
      {
        J(_P_, co, _C_, co) -= intensity * m_imDensityCo.deriv (co, _C_, co);
        J(_P_, co, _T_, co) -= intensity * m_imDensityCo.deriv (co, _T_, co);
      }
    
      if (!m_BoussinesqEnergy)
      {
        if (!m_imDensityCo.constant())
        {
          J(_T_, co, _C_, co)
            -= intensity * m_imDensityCo.deriv(co, _C_, co) * m_imHeatCapacityFluid * u(_T_, co);
          J(_T_, co, _T_, co)
            -= intensity * m_imHeatCapacityFluid * (m_imDensityCo.deriv(co, _T_, co) * u(_T_, co) + m_imDensityCo[co]);
        }
        else
          J(_T_, co, _T_, co) -= intensity * m_imHeatCapacityFluid * m_imDensityCo[co];
      }
      else
        J(_T_, co, _T_, co) -= intensity * m_imHeatCapacityFluid * m_BoussinesqDensity;
    
      if(m_bVolStabDataActive)
        m_spVolStabData->stiff(pElem->vertex(co))
          -= m_BoussinesqFlow? intensity : intensity * m_imDensityCo[co];
    }
  }
  else
  { // heat source
    if (temperature < 0)
      UG_THROW ("FractTHF_FV1: Heat sinks are not implemented.");
    // otherwise nothing to assemble in the matrix
  }
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::ass_dM_elem
(
  LocalVector & d, 
  const LocalVector & u, 
  GridObject * elem, 
  const position_type vCornerCoords [] 
)
{
  TElem * pElem = static_cast<TElem*> (elem);
 
//  call the specific routine:
  if (! m_isFracture)
    this->template bulk_ass_dM_elem<TElem> (d, u, pElem, vCornerCoords);
  else
    this->template fract_ass_dM_elem<TElem> (d, u, pElem, vCornerCoords);
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::bulk_ass_dM_elem
(
  LocalVector & d, 
  const LocalVector & u, 
  TElem * pElem, 
  const position_type vCornerCoords [] 
)
{
  typedef FV1Geometry<TElem, dim> TBulkFVGeom;
  
//  Get finite volume geometry
  const TBulkFVGeom& bulk_geo = GeomProvider<TBulkFVGeom>::get();
  
//  The porosity is constant per element
  number porosity = m_imPorosity [0];
 
//  Loop Sub Control Volumes (SCV)
  for (size_t ip = 0; ip < bulk_geo.num_scv(); ip++)
  {
  //  Get current SCV
    const typename TBulkFVGeom::SCV& scv = bulk_geo.scv(ip);
 
  //  Get associated node
    const int co = scv.node_id();
 
  //  Add to local matrix
    if(m_BoussinesqTransport)
      d(_C_,co) += porosity * u(_C_,co) * scv.volume();
    else
      d(_C_,co) += porosity * m_imDensityCo[ip] * u(_C_,co) * scv.volume();
 
    if(m_BoussinesqFlow)
      d(_P_,co) += porosity * scv.volume();
    else
    {
      d(_P_,co) += porosity * m_imDensityCo[ip] * scv.volume();
      
      if (m_bVolStabDataActive)
      {
      //  Volume stabilization: Add the continuity equation to the transport equation
      //  REMARK: We assume the implicit Euler scheme
        m_spVolStabData->mass(pElem->vertex(co))
          -= porosity * (m_imDensityCo[ip] - m_imOldDensityCo[ip])
                      * scv.volume();
      }
    }
    
    if(m_BoussinesqEnergy)
      d(_T_,co) +=
        (porosity * m_imHeatCapacityFluid * m_BoussinesqDensity
          + (1 - porosity) * m_imHeatCapacitySolid * m_imMassDensitySolid)
            * u(_T_, co) * scv.volume();
    else
      d(_T_,co) +=
        (porosity * m_imHeatCapacityFluid * m_imDensityCo[ip]
          + (1 - porosity) * m_imHeatCapacitySolid * m_imMassDensitySolid)
            * u(_T_, co) * scv.volume();
  }
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::fract_ass_dM_elem
(
  LocalVector & d, 
  const LocalVector & u, 
  TElem * pElem, 
  const position_type vCornerCoords [] 
)
{
//  The porosity (constant per element)
  number porosity = m_imPorosity[0];
  
//  Loop Sub Control Volumes (SCV)
  for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
  {
  //  Get current SCV
    const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv(ip);
 
  //  Fracture width ('/ 2' because every fracture is a boundary of 2 elements)
    number vol = scv.volume() * m_imAperture[0] / 2;
    
  //  Get associated node of the element (not side!)
    const int co = m_innerSideCo [scv.node_id()];
 
  //  Add to local matrix
    if(m_BoussinesqTransport)
      d(_C_,co) += porosity * u(_C_,co) * vol;
    else
      d(_C_,co) += porosity * m_imDensityCo[co] * u(_C_,co) * vol;
 
    if(m_BoussinesqFlow)
      d(_P_,co) += porosity * vol;
    else
    {
      d(_P_,co) += porosity * m_imDensityCo[co] * vol;
      
      if (m_bVolStabDataActive)
      {
      //  Volume stabilization: Add the continuity equation to the transport equation
      //  REMARK: We assume the implicit Euler scheme
        m_spVolStabData->mass(pElem->vertex(co))
          -= porosity * (m_imDensityCo[co] - m_imOldDensityCo[co]) * vol;
      }
    }
    
    if(m_BoussinesqEnergy)
      d(_T_,co) +=
        (porosity * m_imHeatCapacityFluid * m_BoussinesqDensity
          + (1 - porosity) * m_imHeatCapacitySolid * m_imMassDensitySolid)
            * u(_T_, co) * vol;
    else
      d(_T_,co) +=
        (porosity * m_imHeatCapacityFluid * m_imDensityCo[co]
          + (1 - porosity) * m_imHeatCapacitySolid * m_imMassDensitySolid)
            * u(_T_, co) * vol;
  }
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::ass_JM_elem
(
  LocalMatrix & J, 
  const LocalVector & u, 
  GridObject * elem, 
  const position_type vCornerCoords [] 
)
{
  TElem * pElem = static_cast<TElem*> (elem);
  
//  call the specific routine:
  if (! m_isFracture)
    this->template bulk_ass_JM_elem<TElem> (J, u, pElem, vCornerCoords);
  else
    this->template fract_ass_JM_elem<TElem> (J, u, pElem, vCornerCoords);
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::bulk_ass_JM_elem
(
  LocalMatrix & J, 
  const LocalVector & u, 
  TElem * pElem, 
  const position_type vCornerCoords [] 
)
{
  typedef FV1Geometry<TElem, dim> TBulkFVGeom;
  
//  get finite volume geometry
  const TBulkFVGeom& bulk_geo = GeomProvider<TBulkFVGeom>::get();
 
//  The porosity is constant per element
  number porosity = m_imPorosity [0];
 
//  loop Sub Control Volumes (SCV)
  for (size_t ip = 0; ip < bulk_geo.num_scv(); ip++)
  {
  //  get current SCV
    const typename TBulkFVGeom::SCV& scv = bulk_geo.scv(ip);
 
  //  get associated node
    const int co = scv.node_id();
 
  //  Add to local matrix
    if(m_BoussinesqTransport)
      J(_C_, co, _C_, co) += porosity * scv.volume();
    else
      J(_C_, co, _C_, co) +=
        porosity * scv.volume() *
          (m_imDensityCo[ip] + m_imDensityCo.deriv(ip, _C_, co) * u(_C_, co));
 
    if(!m_BoussinesqFlow)
    {
      J(_P_, co, _C_, co) += porosity * m_imDensityCo.deriv(ip, _C_, co) * scv.volume();
      
      if (m_bVolStabDataActive)
      {
      //  Volume stabilization: Add the continuity equation to the transport equation
      //  REMARK: We assume the implicit Euler scheme
        m_spVolStabData->mass(pElem->vertex(co))
          -= porosity * scv.volume() * (m_imDensityCo[ip] - m_imOldDensityCo[ip]);
        
        /*TODO: Do we need this derivative? What we do is subtracting the defect of the flow
         * equation multiplied by _C_ from the transport equation. We differentiate the transport
         * equation, but not this artificial sink. Thus, we should skip this derivative. Otherwise
         * the contribution of the sink will not converge to 0 as soon as we approach the solution
         * of the system. However, in ug3, this derivative is present.
         *
        if(!m_imDensityScv.constant())
          m_spVolStabData->mass(pElem->vertex(co)) -=
            porosity * scv.volume() * m_imDensityCo.deriv(ip, _C_, co)*u(_C_, co);
         */
      }
    }
    //else
      //J(_P_, co, _C_, co) += 0;
 
    if(m_BoussinesqEnergy)
      J(_T_, co, _T_, co) +=
        (porosity * m_imHeatCapacityFluid * m_BoussinesqDensity
        + (1 - porosity) * m_imHeatCapacitySolid * m_imMassDensitySolid)
          * scv.volume();
    else
    {
      J(_T_, co, _T_, co) +=
        (
          (porosity * m_imHeatCapacityFluid * m_imDensityCo[ip]
          + (1 - porosity) * m_imHeatCapacitySolid * m_imMassDensitySolid)
          + porosity * m_imHeatCapacityFluid * (m_imDensityCo.constant() ? 0.0 : m_imDensityCo.deriv(ip, _T_, co))
            * u(_T_, co)
        ) * scv.volume();
 
      J(_T_, co, _C_, co) +=
        porosity * m_imHeatCapacityFluid * (m_imDensityCo.constant() ? 0.0 : m_imDensityCo.deriv(ip, _C_, co))
            * u(_T_, co) * scv.volume();
    }
  }
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::fract_ass_JM_elem
(
  LocalMatrix & J, 
  const LocalVector & u, 
  TElem * pElem, 
  const position_type vCornerCoords [] 
)
{
//  The porosity (constant per element)
  number porosity = m_imPorosity[0];
 
//  loop Sub Control Volumes (SCV)
  for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
  {
  //  get current SCV
    const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv(ip);
 
  //  Fracture width ('/ 2' because every fracture is a boundary of 2 elements)
    number vol = scv.volume() * m_imAperture[0] / 2;
    
  //  Get associated node of the element (not side!)
    const int co = m_innerSideCo [scv.node_id()];
 
  //  Add to local matrix
    if(m_BoussinesqTransport)
      J(_C_, co, _C_, co) += porosity * vol;
    else
    {
      J(_C_, co, _C_, co) +=
        porosity * vol *
          (m_imDensityCo[ip]
           + (m_imDensityCo.constant() ? 0.0 : m_imDensityCo.deriv(co, _C_, co)) * u(_C_, co));
      J(_C_, co, _T_, co) +=
        porosity * vol *
             (m_imDensityCo.constant() ? 0.0 : m_imDensityCo.deriv(co, _T_, co)) * u(_C_, co);
    }
 
    if(!m_BoussinesqFlow){
      J(_P_, co, _C_, co) += porosity *  (m_imDensityCo.constant() ? 0.0 : m_imDensityCo.deriv(co, _C_, co)) * vol;
      J(_P_, co, _T_, co) += porosity *  (m_imDensityCo.constant() ? 0.0 : m_imDensityCo.deriv(co, _T_, co)) * vol;
      
      if (m_bVolStabDataActive)
      {
      //  Volume stabilization: Add the continuity equation to the transport equation
      //  REMARK: We assume the implicit Euler scheme
        m_spVolStabData->mass(pElem->vertex(co))
          -= porosity * vol * (m_imDensityCo[co] - m_imOldDensityCo[co]);
        
        /*TODO: Do we need this derivative? What we do is subtracting the defect of the flow
         * equation multiplied by _C_ from the transport equation. We differentiate the transport
         * equation, but not this artificial sink. Thus, we should skip this derivative. Otherwise
         * the contribution of the sink will not converge to 0 as soon as we approach the solution
         * of the system. However, in ug3, this derivative is present.
         *
        if(!m_imDensityScv.constant())
          m_spVolStabData->mass(pElem->vertex(co)) -=
            porosity * vol * m_imDensityCo.deriv(co, _C_, co)*u(_C_, co);
         */
      }
    }
    //else
      //J(_P_, co, _C_, co) += 0;
      //J(_P_, co, _T_, co) += 0;
 
 
    if(m_BoussinesqEnergy)
      J(_T_,co, _T_, co) += vol *
          (porosity * m_imHeatCapacityFluid * m_BoussinesqDensity
          + (1 - porosity) * m_imHeatCapacitySolid * m_imMassDensitySolid);
    else
    {
      J(_T_, co, _T_, co) += vol *
        (
          (porosity * m_imHeatCapacityFluid * m_imDensityCo[co]
          + (1 - porosity) * m_imHeatCapacitySolid * m_imMassDensitySolid)
          + porosity * m_imHeatCapacityFluid * (m_imDensityCo.constant() ? 0.0 : m_imDensityCo.deriv(co, _T_, co))
            * u(_T_, co)
        );
 
      J(_T_, co, _C_, co) += vol *
            porosity * m_imHeatCapacityFluid * (m_imDensityCo.constant() ? 0.0 : m_imDensityCo.deriv(co, _C_, co))
              * u(_T_, co);
    }
 
    // Remark: Other addings are zero
    //J(_C_, co, _P_, co) += 0;
    //J(_P_, co, _P_, co) += 0;
    //J(_T_, co, _P_, co) += 0;
  }
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::ass_rhs_elem
(
  LocalVector & b, 
  GridObject * elem, 
  const position_type vCornerCoords [] 
)
{
}
 
//  implementation of the export parameters
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::
ex_brine(number vValue[],
         const MathVector<dim> vGlobIP[],
         number time, int si,
         const LocalVector& u,
         GridObject* elem,
         const MathVector<dim> vCornerCoords[],
         const MathVector<dim> vLocIP[],
         const size_t nip,
         bool bDeriv,
         std::vector<std::vector<number> > vvvDeriv[])
{
//  reference element
  typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
 
//  number of shape functions
  static const size_t numSH = ref_elem_type::numCorners;
 
//  get trial space
  LagrangeP1<ref_elem_type> rTrialSpace = Provider<LagrangeP1<ref_elem_type> >::get();
 
//  storage for shape function at ip
  number vShape[numSH];
 
//  loop ips
  for(size_t ip = 0; ip < nip; ++ip)
  {
  //  evaluate at shapes at ip
    rTrialSpace.shapes(vShape, vLocIP[ip]);
 
  //  compute concentration at ip
    vValue[ip] = 0.0;
    for(size_t sh = 0; sh < numSH; ++sh)
      vValue[ip] += u(_C_, sh) * vShape[sh];
 
  //  compute derivative w.r.t. to unknowns iff needed
    if(bDeriv)
      for(size_t sh = 0; sh < numSH; ++sh)
      {
        vvvDeriv[ip][_C_][sh] = vShape[sh];
        vvvDeriv[ip][_P_][sh] = 0.0;
        vvvDeriv[ip][_T_][sh] = 0.0;
      }
  }
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::
ex_pressure(number vValue[],
         const MathVector<dim> vGlobIP[],
         number time, int si,
         const LocalVector& u,
         GridObject* elem,
         const MathVector<dim> vCornerCoords[],
         const MathVector<dim> vLocIP[],
         const size_t nip,
         bool bDeriv,
         std::vector<std::vector<number> > vvvDeriv[])
{
//  reference element
  typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
 
//  number of shape functions
  static const size_t numSH = ref_elem_type::numCorners;
 
//  get trial space
  LagrangeP1<ref_elem_type> rTrialSpace = Provider<LagrangeP1<ref_elem_type> >::get();
 
//  storage for shape function at ip
  number vShape[numSH];
 
//  loop ips
  for(size_t ip = 0; ip < nip; ++ip)
  {
  //  evaluate at shapes at ip
    rTrialSpace.shapes(vShape, vLocIP[ip]);
 
  //  compute concentration at ip
    vValue[ip] = 0.0;
    for(size_t sh = 0; sh < numSH; ++sh)
      vValue[ip] += u(_P_, sh) * vShape[sh];
 
  //  compute derivative w.r.t. to unknowns iff needed
    if(bDeriv)
      for(size_t sh = 0; sh < numSH; ++sh)
      {
        vvvDeriv[ip][_P_][sh] = vShape[sh];
        vvvDeriv[ip][_C_][sh] = 0.0;
        vvvDeriv[ip][_T_][sh] = 0.0;
      }
  }
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::
ex_temperature(number vValue[],
         const MathVector<dim> vGlobIP[],
         number time, int si,
         const LocalVector& u,
         GridObject* elem,
         const MathVector<dim> vCornerCoords[],
         const MathVector<dim> vLocIP[],
         const size_t nip,
         bool bDeriv,
         std::vector<std::vector<number> > vvvDeriv[])
{
//  reference element
  typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
 
//  number of shape functions
  static const size_t numSH = ref_elem_type::numCorners;
 
//  get trial space
  LagrangeP1<ref_elem_type> rTrialSpace = Provider<LagrangeP1<ref_elem_type> >::get();
 
//  storage for shape function at ip
  number vShape[numSH];
 
//  loop ips
  for(size_t ip = 0; ip < nip; ++ip)
  {
  //  evaluate at shapes at ip
    rTrialSpace.shapes(vShape, vLocIP[ip]);
 
  //  compute concentration at ip
    vValue[ip] = 0.0;
    for(size_t sh = 0; sh < numSH; ++sh)
      vValue[ip] += u(_T_, sh) * vShape[sh];
 
  //  compute derivative w.r.t. to unknowns iff needed
    if(bDeriv)
      for(size_t sh = 0; sh < numSH; ++sh)
      {
        vvvDeriv[ip][_T_][sh] = vShape[sh];
        vvvDeriv[ip][_C_][sh] = 0.0;
        vvvDeriv[ip][_P_][sh] = 0.0;
      }
  }
}
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::
ex_darcy_fract
(
  MathVector<dim> vValue[],
  const MathVector<dim> vGlobIP[],
  number time, int si,
  const LocalVector& u,
  GridObject* elem,
  const MathVector<dim> vCornerCoords[],
  const MathVector<dim> vLocIP[],
  const size_t nip,
  bool bDeriv,
  std::vector<std::vector<MathVector<dim> > > vvvDeriv[]
)
{
  if (m_isFracture)
  {
    typedef StdLinConsistentGravity<low_dim> TConsGravity;
 
  //  this evaluation is implemented only for the standard set of the IPs
    if (nip != m_pFractGeo->num_scvf())
      UG_THROW ("FractTHF_FV1: The Darcy velocity export parameter is only implemented for the standard set of IPs.");
 
  //  number of shape functions
      const size_t numSh = m_pFractGeo->num_sh();
 
    MathVector<dim> Vel;
    MathVector<dim> Vel_c[maxFractSideCorners], Vel_p[maxFractSideCorners], Vel_T[maxFractSideCorners];
    number vViscosity_c [maxFractSideCorners], * pViscosity_c;
    number vViscosity_T [maxFractSideCorners], * pViscosity_T;
 
  //  Get the corner values
    number vDensity [maxFractSideCorners];
    number vPressure [maxFractSideCorners];
    for (size_t sh = 0; sh < numSh; sh++)
    {
      size_t co = m_innerSideCo[sh];
      vDensity[sh] = m_imDensityCo[co];
      vPressure[sh] = u(_P_, co);
    }
 
  //  Consistent gravity and its derivative at corners
    TConsGravity ConsGravityMethod;
      MathVector<low_dim> vConsGravity [maxFractSideCorners];
      MathVector<low_dim> vConsGravity_c [maxFractSideCorners][maxFractSideCorners];
      MathVector<low_dim> vConsGravity_T [maxFractSideCorners][maxFractSideCorners];
 
  //  Prepare the consistent gravity
    try
    {
      ConsGravityMethod.template prepare<dim>
        (vConsGravity, numSh, m_pFractGeo->corners(), vDensity, m_Gravity);
 
        if(bDeriv)
        {
          number co_density [maxFractSideCorners];
          memset (co_density, 0, numSh * sizeof (number));
          for (size_t sh = 0; sh < numSh; sh++)
          { //TODO: use scv, not merely shape id here!
            size_t co = m_innerSideCo[sh];
          if(!m_imDensityCo.constant())
            co_density[sh] = m_imDensityCo.deriv (co, _C_, co); // we assume the completely local dependence of \f$\rho\f$ of \f$c\f$
            ConsGravityMethod.template prepare<dim>
              (vConsGravity_c[sh], numSh, m_pFractGeo->corners(), co_density, m_Gravity);
            co_density[sh] = 0.0;
          }
          for (size_t sh = 0; sh < numSh; sh++)
          { //TODO: use scv, not merely shape id here!
            size_t co = m_innerSideCo[sh];
          if(!m_imDensityCo.constant())
            co_density[sh] = m_imDensityCo.deriv (co, _T_, co); // we assume the completely local dependence of \f$\rho\f$ of \f$c\f$
            ConsGravityMethod.template prepare<dim>
              (vConsGravity_T[sh], numSh, m_pFractGeo->corners(), co_density, m_Gravity);
            co_density[sh] = 0.0;
          }
        }
    }
    UG_CATCH_THROW ("FractTHF_FV1::fract_ass_dA_elem: Cannot prepare Consistent Gravity.");
 
  //  Loop the ips
    for (size_t ip = 0; ip < nip; ip++)
    {
    //  Darcy velocity to be filled
      MathVector<dim>& DarcyVel = vValue[ip];
 
    //  Compute the Darcy velocity at the ips and the convection shapes
        this->template compute_ip_Darcy_velocity<TFractFVGeom, TConsGravity>
            (Vel, ip, *m_pFractGeo, ConsGravityMethod, vConsGravity, vPressure,
              m_imViscosityIP[ip]);
        VecScale (DarcyVel, Vel, m_imFractPermeability[0]);
        
    //  Compute the derivatives
        if(bDeriv)
        {
      //  Derivatives to be filled
        MathVector<dim>* DarcyVel_c = &vvvDeriv[ip][_C_][0];
        MathVector<dim>* DarcyVel_p = &vvvDeriv[ip][_P_][0];
        MathVector<dim>* DarcyVel_T = &vvvDeriv[ip][_T_][0];
        
      //  Clear all the derivatives (as we fill up not all of them):
        for (size_t co = 0; co < vvvDeriv[ip][_C_].size(); co++)
          DarcyVel_c[co] = 0.0;
        for (size_t co = 0; co < vvvDeriv[ip][_P_].size(); co++)
          DarcyVel_p[co] = 0.0;
        for (size_t co = 0; co < vvvDeriv[ip][_T_].size(); co++)
          DarcyVel_T[co] = 0.0;
        
      //  Compute the derivative of the viscosity
        if (m_imViscosityIP.constant ())
        {
          pViscosity_c = NULL;
          pViscosity_T = NULL;
        }
        else
        {
          for (size_t sh = 0; sh < numSh; sh++)
          {
              size_t co = m_innerSideCo[sh];
            vViscosity_c[sh] = m_imViscosityIP.deriv(ip, _C_, co);
            vViscosity_T[sh] = m_imViscosityIP.deriv(ip, _T_, co);
          }
          pViscosity_c = vViscosity_c;
          pViscosity_T = vViscosity_T;
        }
 
      // Compute derivatives of the Darcy velocity
        this->template compute_J_ip_Darcy_velocity <TFractFVGeom, TConsGravity, maxFractSideCorners>
                (Vel, Vel_c, Vel_p, Vel_T, ip, *m_pFractGeo, ConsGravityMethod,
                  vConsGravity_c, vConsGravity_T, vPressure,
                    m_imViscosityIP[ip], pViscosity_c, pViscosity_T);
        for (size_t sh = 0; sh < numSh; sh++)
        {
            size_t co = m_innerSideCo[sh];
          VecScale (DarcyVel_c[co], Vel_c[sh], m_imFractPermeability[0]);
            VecScale (DarcyVel_p[co], Vel_p[sh], m_imFractPermeability[0]);
            VecScale (DarcyVel_T[co], Vel_T[sh], m_imFractPermeability[0]);
        }
        }
    }
  }
  else
    UG_THROW ("FractTHF_FV1: The Darcy velocity export parameter is currently implemented only for fractures.");
}
 
 
template<typename TDomain>
template<typename TElem>
void FractTHF_FV1<TDomain>::
ex_darcy_ortho_fract
(
  number vValue[],
  const MathVector<dim> vGlobIP[],
  number time, int si,
  const LocalVector& u,
  GridObject* elem,
  const MathVector<dim> vCornerCoords[],
  const MathVector<dim> vLocIP[],
  const size_t nip,
  bool bDeriv,
  std::vector<std::vector<number> > vvvDeriv[]
)
{
  if (m_isFracture)
  {
  //  this evaluation is implemented only for the standard set of the IPs
    if (nip != m_pFractGeo->num_scv())
      UG_THROW ("FractDDF_FV1: The Darcy velocity export parameter is only implemented for the standard set of IPs.");
 
    const number half_fr_width = m_imAperture[0] / 2;
  
  //  loop over the corners of the inner side
    for (size_t ip = 0; ip < nip; ip++)
    {
    //  Get current SCV
      const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv(ip);
      number& orthVelocity = vValue[ip];
 
    //  Get associated node of the element (not side!)
      const int co = m_innerSideCo [scv.node_id()];
    
    //  Get the corner values
      const number orthP_f = u(_P_, co);
      const number orthP_m = u(_P_, m_assCo[co]);
      const number fractDensity = m_imDensityCo[co];
      const number orthDensity = m_imDensityCo [m_assCo[co]];
      const number orthViscosity = m_imViscosityCo [m_assCo[co]];
  
    //  Flux velocity between the fracture and the bulk medium
      orthVelocity = (m_orthGravity * (orthDensity - fractDensity)
        - (orthP_m - orthP_f) / half_fr_width)
          * m_imOrthoPermeability[0];
        
    //  Compute the derivatives
        if(bDeriv)
        {
      //  Derivatives to be filled
        number* DarcyVel_c = &vvvDeriv[ip][_C_][0]; // array indexed by the shape
        number* DarcyVel_p = &vvvDeriv[ip][_P_][0]; // array indexed by the shape
        number* DarcyVel_T = &vvvDeriv[ip][_T_][0]; // array indexed by the shape
        
      //  Clear all the derivatives (as we fill up not all of them):
        for (size_t i = 0; i < vvvDeriv[ip][_C_].size(); i++) DarcyVel_c[i] = 0;
        for (size_t i = 0; i < vvvDeriv[ip][_P_].size(); i++) DarcyVel_p[i] = 0;
        for (size_t i = 0; i < vvvDeriv[ip][_T_].size(); i++) DarcyVel_p[i] = 0;
      
      //  Compute the derivatives:
        number D_fractDensity_c, D_orthDensity_c, D_fractDensity_T, D_orthDensity_T;
        if (! m_imDensityCo.constant())
        {
          D_fractDensity_c = m_imDensityCo.deriv (co, _C_, co);
          D_fractDensity_T = m_imDensityCo.deriv (co, _T_, co);
          D_orthDensity_c = m_imDensityCo.deriv (m_assCo[co], _C_, m_assCo[co]);
          D_orthDensity_T = m_imDensityCo.deriv (m_assCo[co], _T_, m_assCo[co]);
        }
        else
          D_fractDensity_c = D_fractDensity_T = D_orthDensity_c = D_orthDensity_T = 0;
  
        const number orthViscosity = m_imViscosityCo [m_assCo[co]];
        number D_orthViscosity_c, D_orthViscosity_T;
        if (! m_imViscosityCo.constant())
        {
          D_orthViscosity_c = m_imViscosityCo.deriv (m_assCo[co], _C_, m_assCo[co]);
          D_orthViscosity_T = m_imViscosityCo.deriv (m_assCo[co], _T_, m_assCo[co]);
        }
        else
          D_orthViscosity_c = D_orthViscosity_T = 0;
 
        number D_orthVelocity [3], D_orthVelocity_fr [3]; // 3 components: w.r.t. _C_, _P_ and _T_ resp.
    
        D_orthVelocity [_C_] = (m_orthGravity * D_orthDensity_c
          * m_imOrthoPermeability[0] * orthViscosity
            - orthVelocity * D_orthViscosity_c) / orthViscosity / orthViscosity;
    
        D_orthVelocity_fr [_C_] = - m_orthGravity * D_fractDensity_c
          * m_imOrthoPermeability[0] / orthViscosity;
    
        D_orthVelocity [_P_] = - m_imOrthoPermeability[0]
          / half_fr_width / orthViscosity;
    
        D_orthVelocity_fr [_P_] = - D_orthVelocity [_P_];
    
        D_orthVelocity [_T_] = (m_orthGravity * D_orthDensity_T
          * m_imOrthoPermeability[0] * orthViscosity
            - orthVelocity * D_orthViscosity_T) / orthViscosity / orthViscosity;
    
        D_orthVelocity_fr [_T_] = - m_orthGravity * D_fractDensity_T
          * m_imOrthoPermeability[0] / orthViscosity;
        
        DarcyVel_c[co] = D_orthVelocity_fr[_C_];
        DarcyVel_p[co] = D_orthVelocity_fr[_P_];
        DarcyVel_T[co] = D_orthVelocity_fr[_T_];
        
        DarcyVel_c[m_assCo[co]] = D_orthVelocity[_C_];
        DarcyVel_p[m_assCo[co]] = D_orthVelocity[_P_];
        DarcyVel_T[m_assCo[co]] = D_orthVelocity[_T_];
        }
        
    //  The viscosity:
          orthVelocity /= orthViscosity;
    }
  }
  else
    UG_THROW ("FractTHF_FV1: The orthogonal Darcy velocity export parameter is implemented only for fractures.");
}
 
//  register assembling functions
 
template<typename TDomain>
template<typename TElem> // the element to register for
void FractTHF_FV1<TDomain>::register_loc_discr_func ()
{
  static const ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
  static const int refDim = reference_element_traits<TElem>::dim;
  
  this->clear_add_fct(id);
  
  this->set_prep_elem_loop_fct(id, & this_type::template prepare_element_loop<TElem>);
  this->set_prep_elem_fct   (id, & this_type::template prepare_element<TElem>);
  this->set_fsh_elem_loop_fct (id, & this_type::template finish_element_loop<TElem>);
  this->set_add_jac_A_elem_fct(id, & this_type::template ass_JA_elem<TElem>);
  this->set_add_jac_M_elem_fct(id, & this_type::template ass_JM_elem<TElem>);
  this->set_add_def_A_elem_fct(id, & this_type::template ass_dA_elem<TElem>);
  this->set_add_def_M_elem_fct(id, & this_type::template ass_dM_elem<TElem>);
  this->set_add_rhs_elem_fct  (id, & this_type::template ass_rhs_elem<TElem>);
  
  m_exBrine->template set_fct<this_type, refDim> (id, this, &this_type::template ex_brine<TElem>);
  m_exPressure->template set_fct<this_type, refDim> (id, this, &this_type::template ex_pressure<TElem>);
  m_exTemperature->template set_fct<this_type, refDim> (id, this, &this_type::template ex_temperature<TElem>);
  m_exFractDarcyVel->template set_fct<this_type,refDim>(id, this, &this_type::template ex_darcy_fract<TElem>);
  m_exOrthoFractDarcyVel->template set_fct<this_type,refDim>(id, this, &this_type::template ex_darcy_ortho_fract<TElem>);
}
 
template <typename TDomain>
FractTHF_FV1<TDomain>::FractTHF_FV1
(
  const char* functions, 
  const char* subsets 
)
: IElemDisc<TDomain> (functions, subsets)
{
  init();
};
 
template <typename TDomain>
FractTHF_FV1<TDomain>::FractTHF_FV1
(
    const std::vector<std::string>& vFct,
    const std::vector<std::string>& vSubset
)
: IElemDisc<TDomain> (vFct, vSubset)
{
  init();
};
 
template <typename TDomain>
void FractTHF_FV1<TDomain>::init()
{
  set_boussinesq (false);
  
  m_imDensityCo.set_comp_lin_defect (false);
  m_imOldDensityCo.set_comp_lin_defect (false);
  m_imDensityIP.set_comp_lin_defect (false),
  m_imViscosityCo.set_comp_lin_defect (false);
  m_imViscosityIP.set_comp_lin_defect (false);
  m_imConstGravity.set_comp_lin_defect (false);
  m_imPorosity.set_comp_lin_defect (false);
  m_imPermeability .set_comp_lin_defect(false);
  m_imDiffusion.set_comp_lin_defect (false);
  m_imThermalConductivity.set_comp_lin_defect (false);
  m_imAperture.set_comp_lin_defect (false);
  m_imFractPermeability.set_comp_lin_defect (false);
  m_imOrthoPermeability.set_comp_lin_defect (false);
  m_imFractDiffusion.set_comp_lin_defect (false);
  m_imOrthoDiffusion.set_comp_lin_defect (false);
  m_imFractThermCond.set_comp_lin_defect (false);
  m_imOrthoThermCond.set_comp_lin_defect (false);
  m_spUpwind = make_sp(new ConvectionShapesFullUpwind<dim>);
  m_spUpwindEnergy = make_sp(new ConvectionShapesFullUpwind<dim>);
  m_imFractDarcyVelIP.set_comp_lin_defect (false);
 
  m_spVolStabData = SPNULL;
  
  std::string functions;
  for(size_t i = 0; i < this->symb_fcts().size(); ++i)
  {
    if(i > 0) functions.append(",");
    functions.append(this->symb_fcts()[i]);
  }
 
  m_exBrine = make_sp (new DataExport<number, dim> (functions.c_str()));
  m_exPressure = make_sp (new DataExport<number, dim> (functions.c_str()));
  m_exTemperature = make_sp (new DataExport<number, dim> (functions.c_str()));
  m_exFractDarcyVel = make_sp (new DataExport<MathVector<dim>, dim>(functions.c_str()));
  m_exOrthoFractDarcyVel = make_sp (new DataExport<number, dim>(functions.c_str()));
 
//  check number of functions
  if (this->num_fct () != 3)
    UG_THROW ("Wrong number of functions: The ElemDisc 'FractTHF_FV1'"
          " needs exactly 3 symbolic function"
          " (one for the mass fraction and one for the pressure).");
 
//  register imports
  
  this->register_import (m_imDensityCo);
  this->register_import (m_imOldDensityCo);
  this->register_import (m_imDensityIP);
  this->register_import (m_imViscosityCo);
  this->register_import (m_imViscosityIP);
  this->register_import (m_imConstGravity);
  
  this->register_import (m_imPorosity);
  
  this->register_import (m_imPermeability);
  this->register_import (m_imDiffusion);
  this->register_import (m_imThermalConductivity);
  
  this->register_import (m_imAperture);
  this->register_import (m_imFractPermeability);
  this->register_import (m_imOrthoPermeability);
  this->register_import (m_imFractDiffusion);
  this->register_import (m_imOrthoDiffusion);
  this->register_import (m_imFractThermCond);
  this->register_import (m_imOrthoThermCond);
  
  this->register_import(m_imFractDarcyVelIP);
 
//  connect to own export
  m_imFractDarcyVelIP.set_data(m_exFractDarcyVel);
 
//  register the local assembling functions
  boost::mpl::for_each<AssembleElemList> (RegisterLocalDiscr (this));
 
}
 
} // namespace d3f
} // end namespace ug
 
/* End of File */