fract_pf_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
 * 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_pf_fv1_impl.h - class member implementation of the discretization of
 * the FV-discretization of the pressure-driven flow in fractured porous media.
 */
/* 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 FractPF_FV1<TDomain>::prepare_setting
(
  const std::vector<LFEID> & vLfeID,
  bool bNonRegular
)
{
//  check the grid
  if (bNonRegular)
    UG_THROW ("ERROR in FractPF_FV1:"
        " The discretization does not support hanging nodes.\n");
 
//  check number of the components
  if (vLfeID.size () != 1)
    UG_THROW ("FractPF_FV1: The pressure driven flow requires one coponent 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))
    UG_THROW ("FractPF_FV1: This discretization works with the LagrangeP1-elements only.");
}
 
// assembling
 
template<typename TDomain>
template<typename TElem, typename TBulkFVGeom>
void FractPF_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 ("FractPF_FV1: Fracture manager not closed");
    m_isFracture = m_spFractManager->contains (si);
  }
  
//  get the gravity
  if (!m_imConstGravity.data_given())
    UG_THROW ("FractPF_FV1: Missing Import 'gravity'.");
  if (m_imConstGravity.constant())
    (*m_imConstGravity.user_data()) (m_Gravity, MathVector<dim>(), 0.0, 0);
  else
    UG_THROW ("FractPF_FV1: Gravity must be constant.");
  
//  call the specific routine:
  if (! m_isFracture)
    this->template bulk_prepare_element_loop<TElem, TBulkFVGeom> (roid, si);
  else
    this->template fract_prepare_element_loop<TElem> (roid, si);
}
 
template<typename TDomain>
template<typename TElem, typename TBulkFVGeom>
void FractPF_FV1<TDomain>::bulk_prepare_element_loop
(
  ReferenceObjectID roid, 
  int si 
)
{
//  invalidate the fracture subset imports for this (full-dim.) subset
  m_imAperture.set_data (SPNULL);
  m_imFractPermeability.set_data (SPNULL);
  m_imOrthoPermeability.set_data (SPNULL);
  
//  check the imports
  if (!m_imPermeability.data_given())
    UG_THROW ("FractPF_FV1: Missing Import 'full-dim. permeability'.");
  
//  get the FV geometry
  TBulkFVGeom& bulk_geo = GeomProvider<TBulkFVGeom>::get ();
  static const int refDim = TElem::dim;
  
//  set local IP coordinates
  const MathVector<refDim>* coe_local = bulk_geo.coe_local ();
  m_imPermeability.template set_local_ips<refDim> (coe_local, 1, false);
}
 
template<typename TDomain>
template<typename TElem>
void FractPF_FV1<TDomain>::fract_prepare_element_loop
(
  ReferenceObjectID roid, 
  int si 
)
{
//  initialize the pointer to the FV geometry for fracture elements
  m_pFractGeo = &GeomProvider<TFractFVGeom>::get (LFEID (LFEID::LAGRANGE, low_dim, 1), 1);
  
//  enable the fracture parameters for this subset
  m_imAperture.set_data (m_spApertureData);
  m_imFractPermeability.set_data (m_spFractPermeabilityData);
  m_imOrthoPermeability.set_data (m_spOrthoPermeabilityData);
  
//  check the imports
  if (!m_imAperture.data_given ())
    UG_THROW ("FractPF_FV1: Missing Import 'fracture width (aperture)'.");
  if (!m_imPermeability.data_given ())
    UG_THROW ("FractPF_FV1: Missing Import 'full-dim. permeability' (for the fract.-bulk interface interaction).");
  if (!m_imFractPermeability.data_given ())
    UG_THROW ("FractPF_FV1: Missing Import 'permeability along fracture'.");
  /* If m_imOrthoPermeability is not given, we use m_imPermeability. */
}
 
template<typename TDomain>
template<typename TElem, typename TBulkFVGeom>
void FractPF_FV1<TDomain>::finish_element_loop ()
{
}
 
template<typename TDomain>
template<typename TElem, typename TBulkFVGeom>
void FractPF_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, TBulkFVGeom> (u, pElem, vCornerCoords);
  else
    this->template fract_prepare_element<TElem> (u, pElem, vCornerCoords);
}
 
template<typename TDomain>
template<typename TElem, typename TBulkFVGeom>
void FractPF_FV1<TDomain>::bulk_prepare_element
(
  const LocalVector & u, 
  TElem * elem, 
  const position_type vCornerCoords[] 
)
{
//  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("FractPF_FV1: Cannot update the Finite Volume Geometry for a bulk element.");
 
//  set global positions for user data
  const MathVector<dim>* coe_global = bulk_geo.coe_global ();
  m_imPermeability.set_global_ips (coe_global, 1);
}
 
template<typename TDomain>
template<typename TElem>
void FractPF_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("FractPF_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("FractPF_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("FractPF_FV1: Cannot transform local side coordinates to local element coordinates in a fracture element.");
  
//  get the geometric data of the outer side
  MathVector<dim> loc_center_of_outer_sd, glob_center_of_outer_sd;
  loc_center_of_outer_sd = 0.0; glob_center_of_outer_sd = 0.0;
  for (size_t co = 0; co < m_numFractCo; co++)
  {
    loc_center_of_outer_sd += rRefElem.corner (m_outerSideCo [co]);
    glob_center_of_outer_sd += vCornerCoords [m_outerSideCo [co]];
  }
  loc_center_of_outer_sd /= m_numFractCo;
  glob_center_of_outer_sd /= m_numFractCo;
  
//  set local IP coordinates
  m_imAperture.template set_local_ips<dim> (&m_elem_loc_coe, 1);
  m_imFractPermeability.template set_local_ips<dim> (&m_elem_loc_coe, 1);
  if (m_imOrthoPermeability.data_given ())
    m_imOrthoPermeability.template set_local_ips<dim> (&m_elem_loc_coe, 1);
  else
    m_imPermeability.template set_local_ips<dim> (&m_elem_loc_coe, 1);
  
//  set global positions for user data
  const MathVector<dim>* coe_global = m_pFractGeo->coe_global ();
  m_imAperture.set_global_ips (coe_global, 1);
  m_imFractPermeability.set_global_ips (coe_global, 1);
  if (m_imOrthoPermeability.data_given ())
    m_imOrthoPermeability.set_global_ips (coe_global, 1);
  else
    m_imPermeability.set_global_ips (coe_global, 1);
  
//  set the correct subdomain for the bulk-fracture permeability interface
//TODO: This would be a good idea to get the permeability from the interface
//  if (! m_imPermeability.data()->constant ()) // constant data are the same everywhere!
//    m_imPermeability.data()->set_subset
//      (this->subset_handler().get_subset_index (m_outerFractSide));
//END TODO
  
//  compute the outer normal for the fracture side
  SideNormal<ref_elem_type, dim> (m_unitOuterNormal, m_outerFractSideIdx, vCornerCoords);
  const number outerNormalNorm = VecLength (m_unitOuterNormal);
  if (outerNormalNorm < 1e-32)
    UG_THROW ("FractPF_FV1: Cannot get the normal to a fracture.")
  m_unitOuterNormal /= outerNormalNorm;
}
 
template<typename TDomain>
template<typename TFVGeom>
void FractPF_FV1<TDomain>::compute_ip_velocity
(
  MathVector<dim>& Vel, 
  size_t ip, 
  const TFVGeom& geo, 
    number vPressure[] 
)
{
  const typename TFVGeom::SCVF& scvf = geo.scvf(ip);
 
// The pressure-gradient part:
  Vel = 0.0;
  for (size_t sh = 0; sh < scvf.num_sh(); sh++)
    VecScaleAppend (Vel, - vPressure[sh], scvf.global_grad(sh));
 
// The viscosity factor:
  Vel /= m_viscosity;
}
 
template<typename TDomain>
template<typename TElem, typename TBulkFVGeom>
void FractPF_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, TBulkFVGeom> (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, typename TBulkFVGeom>
void FractPF_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;
  
//  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 ();
 
//  Get the corner pressure
  number vPressure [ref_elem_type::numCorners];
  for (size_t sh = 0; sh < numSh; sh++)
    vPressure [sh] = u (_P_, sh);
  
//  Loop Sub Control Volume Faces (SCVF)
  MathVector<dim> Vel, DarcyVel;
  for (size_t ip = 0; ip < numScvf; ip++)
  {
  //  Compute the Darcy velocity at the ips and the convection shapes
      this->template compute_ip_velocity<TBulkFVGeom> (Vel, ip, bulk_geo, vPressure);
      MatVecMult (DarcyVel, m_imPermeability[0], Vel);
      
  //  Get current SCVF
    const typename TBulkFVGeom::SCVF& scvf = bulk_geo.scvf (ip);
  
  //  Compute flux
    const number flux = VecDot (DarcyVel, scvf.normal ()) * m_density;
 
  //  Add contribution to flow equation
    d (_P_, scvf.from ()) += flux;
    d (_P_, scvf.to ())   -= flux;
  }
}
 
template<typename TDomain>
template<typename TElem>
void FractPF_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;
  
    const size_t numSh = m_pFractGeo->num_sh ();
  const size_t numScvf = m_pFractGeo->num_scvf ();
 
//  Get the corner values
  number vPressure [maxFractSideCorners];
  for (size_t sh = 0; sh < numSh; sh++)
    vPressure [sh] = u (_P_, m_innerSideCo [sh]);
  
//  Loop Sub Control Volume Faces (SCVF)
  MathVector<dim> Vel;
  for (size_t ip = 0; ip < numScvf; ip++)
  {
  //  Compute the Darcy velocity at the ips and the convection shapes
    this->template compute_ip_velocity<TFractFVGeom> (Vel, ip, *m_pFractGeo, vPressure);
    
  //  Get current SCVF
    const typename TFractFVGeom::SCVF& scvf = m_pFractGeo->scvf (ip);
 
  //  Compute flux
    number flux = VecDot (Vel, scvf.normal ()) * m_imFractPermeability[0]
            * m_density * half_fr_width;
 
  //  Add contribution to flow equation
    d (_P_, m_innerSideCo [scvf.from ()]) += flux;
    d (_P_, m_innerSideCo [scvf.to ()])   -= flux;
  }
}
 
template<typename TDomain>
template<typename TElem>
void FractPF_FV1<TDomain>::fract_bulk_ass_dA_elem
(
  LocalVector & d, 
  const LocalVector & u, 
  TElem * pElem, 
  const position_type vCornerCoords[] 
)
{
  const number half_fr_width = m_imAperture[0] / 2;
  
//  the permeability of the fracture-bulk interface
  number orthPermeability;
  if (m_imOrthoPermeability.data_given ())
    orthPermeability = m_imOrthoPermeability[0];
  else
  {
    MathVector<dim> orthPermeabilityVec;
    MatVecMult (orthPermeabilityVec, m_imPermeability[0], m_unitOuterNormal);
    orthPermeability = VecDot (orthPermeabilityVec, m_unitOuterNormal);
  }
  
//  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 orthP_f = u (_P_, co);
    const number orthP_m = u (_P_, m_assCo [co]);
  
  //  Flux velocity between the fracture and the bulk medium
    const number orthVelocity = - (orthP_m - orthP_f) / half_fr_width
                * orthPermeability / m_viscosity;
 
        const number flux = orthVelocity * m_density * scv.volume();
        
        d (_P_, m_assCo [co]) -= flux;
        d (_P_, co)           += flux;
  }
}
 
template<typename TDomain>
template<typename TFVGeom>
void FractPF_FV1<TDomain>::compute_J_ip_velocity
(
  MathVector<dim> Vel_p[], 
  size_t ip, 
  const TFVGeom& geo 
)
{
  const typename TFVGeom::SCVF& scvf = geo.scvf (ip);
    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++)
    VecScale (Vel_p [sh], scvf.global_grad (sh), -1 / m_viscosity);
}
 
template<typename TDomain>
template<typename TElem, typename TBulkFVGeom>
void FractPF_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, TBulkFVGeom> (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, typename TBulkFVGeom>
void FractPF_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;
 
  //  get finite volume geometry
  const TBulkFVGeom& bulk_geo = GeomProvider<TBulkFVGeom>::get();
 
  const size_t numScvf = bulk_geo.num_scvf();
 
//  Loop Sub Control Volume Faces (SCVF)
  MathVector<dim> Vel_p [ref_elem_type::numCorners], DarcyVel_p;
  for (size_t ip = 0; ip < numScvf; ip++)
  {
  //  Get current SCVF
    const typename TBulkFVGeom::SCVF& scvf = bulk_geo.scvf (ip);
 
  // compute derivatives of the Darcy velocity
      this->template compute_J_ip_velocity<TBulkFVGeom> (Vel_p, ip, bulk_geo);
 
  //  add the flux contribution to the stiffness matrix
    for (size_t sh = 0; sh < scvf.num_sh (); sh++)
    {
      MatVecMult (DarcyVel_p, m_imPermeability[0], Vel_p [sh]);
      number DFlux_p = VecDot (DarcyVel_p, scvf.normal ()) * m_density;
      
      J (_P_, scvf.from (), _P_, sh) += DFlux_p;
      J (_P_, scvf.to (),   _P_, sh) -= DFlux_p;
    }
  }
}
 
template<typename TDomain>
template<typename TElem>
void FractPF_FV1<TDomain>::fract_ass_JA_elem
(
  LocalMatrix & J, 
  const LocalVector & u, 
  TElem * pElem, 
  const position_type vCornerCoords[] 
)
{
  const number half_fr_width = m_imAperture[0] / 2;
  
  const size_t numScvf = m_pFractGeo->num_scvf ();
 
//  Loop Sub Control Volume Faces (SCVF)
  MathVector<dim> Vel_p [maxFractSideCorners];
  for (size_t ip = 0; ip < numScvf; ip++)
  {
  //  get current SCVF
    const typename TFractFVGeom::SCVF& scvf = m_pFractGeo->scvf(ip);
 
    //  compute the Jacobian of the Darcy velocity
    this->template compute_J_ip_velocity<TFractFVGeom> (Vel_p, ip, *m_pFractGeo);
  
  //  add the flux contribution to the stiffness matrix
    for (size_t sh = 0; sh < scvf.num_sh(); sh++)
    {
    //  get the derivative of the flux
      number DFlux_p = VecDot (Vel_p [sh], scvf.normal ()) * m_imFractPermeability[0]
                * m_density * half_fr_width;
      
    //  get the connection indices
      size_t co = m_innerSideCo [sh];
      size_t co_from = m_innerSideCo [scvf.from ()];
      size_t co_to = m_innerSideCo [scvf.to ()];
      
    //  add to the matrix
      J (_P_, co_from, _P_, co) += DFlux_p;
      J (_P_, co_to,   _P_, co) -= DFlux_p;
    }
  }
}
 
template<typename TDomain>
template<typename TElem>
void FractPF_FV1<TDomain>::fract_bulk_ass_JA_elem
(
  LocalMatrix & J, 
  const LocalVector & u, 
  TElem * pElem, 
  const position_type vCornerCoords[] 
)
{
  const number half_fr_width = m_imAperture[0] / 2;
  
//  the permeability of the fracture-bulk interface
  number orthPermeability;
  if (m_imOrthoPermeability.data_given ())
    orthPermeability = m_imOrthoPermeability[0];
  else
  {
    MathVector<dim> orthPermeabilityVec;
    MatVecMult (orthPermeabilityVec, m_imPermeability[0], m_unitOuterNormal);
    orthPermeability = VecDot (orthPermeabilityVec, m_unitOuterNormal);
  }
  
//  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 () * m_density;
 
  //  get associated node of the element (not side!)
    const int co = m_innerSideCo [scv.node_id ()];
    
  //  Flux velocity between the fracture and the bulk medium
        const number D_orthFlux_fr = orthPermeability / m_viscosity / half_fr_width * s;
    
    J(_P_, m_assCo [co], _P_, m_assCo [co]) += D_orthFlux_fr;
    J(_P_, co,           _P_, m_assCo [co]) -= D_orthFlux_fr;
    
    J(_P_, m_assCo [co], _P_, co)           -= D_orthFlux_fr;
    J(_P_, co,           _P_, co)           += D_orthFlux_fr;
  }
}
 
template<typename TDomain>
template<typename TElem, typename TBulkFVGeom>
void FractPF_FV1<TDomain>::ass_dM_elem
(
  LocalVector & d, 
  const LocalVector & u, 
  GridObject * elem, 
  const position_type vCornerCoords[] 
)
{
  // Nothing to do: FractPF_FV1 discretizes a stationary equation. No mass matrix!
}
 
template<typename TDomain>
template<typename TElem, typename TBulkFVGeom>
void FractPF_FV1<TDomain>::ass_JM_elem
(
  LocalMatrix & J, 
  const LocalVector & u, 
  GridObject * elem, 
  const position_type vCornerCoords[] 
)
{
  // Nothing to do: FractPF_FV1 discretizes a stationary equation. No mass matrix!
}
 
template<typename TDomain>
template<typename TElem, typename TBulkFVGeom>
void FractPF_FV1<TDomain>::ass_rhs_elem
(
  LocalVector & b, 
  GridObject * elem, 
  const position_type vCornerCoords[] 
)
{
  TElem * pElem = static_cast<TElem*> (elem);
 
//  call the specific routine:
  if (! m_isFracture)
  {
    this->template bulk_ass_rhs_elem<TElem, TBulkFVGeom> (b, pElem, vCornerCoords);
  }
  else
  {
    this->template fract_ass_rhs_elem<TElem> (b, pElem, vCornerCoords);
  }
}
 
template<typename TDomain>
template<typename TElem, typename TBulkFVGeom>
void FractPF_FV1<TDomain>::bulk_ass_rhs_elem
(
  LocalVector & b, 
  TElem * elem, 
  const position_type vCornerCoords[] 
)
{
//  Get finite volume geometry
  const TBulkFVGeom& bulk_geo = GeomProvider<TBulkFVGeom>::get ();
 
  const size_t numScvf = bulk_geo.num_scvf ();
 
//  the gravity part
//  REMARK: Note that for for inner nodes, these contributions shoul sum up
//  to 0 (because the density is constant). The only non-zero gravity
//  right-hand sides may remain at the Neumann (or similar) boundaries:
  MathVector<dim> Vel, DarcyVel;
  for (size_t ip = 0; ip < numScvf; ip++)
  {
  //  Get current SCVF
    const typename TBulkFVGeom::SCVF& scvf = bulk_geo.scvf (ip);
  
  //  Compute K (rho * g) / mu
    VecScale (Vel, m_Gravity, m_density / m_viscosity);
      MatVecMult (DarcyVel, m_imPermeability[0], Vel);
 
  //  Compute flux
    number flux = VecDot (DarcyVel, scvf.normal ()) * m_density;
 
  //  Add contribution to the rhs of the flow equation
    b (_P_, scvf.from ()) -= flux;
    b (_P_, scvf.to ())   += flux;
  }
  
//  mass sources
}
 
template<typename TDomain>
template<typename TElem>
void FractPF_FV1<TDomain>::fract_ass_rhs_elem
(
  LocalVector & b, 
  TElem * elem, 
  const position_type vCornerCoords[] 
)
{
  const number half_fr_width = m_imAperture[0] / 2;
  
  const size_t numScvf = m_pFractGeo->num_scvf ();
 
//  the gravity part
//  REMARK: Note that for for inner nodes, these contributions shoul sum up
//  to 0 (because the density is constant). The only non-zero gravity
//  right-hand sides may remain at the Neumann (or similar) boundaries:
  MathVector<dim> DarcyVel;
  for (size_t ip = 0; ip < numScvf; ip++)
  {
  //  Get current SCVF
    const typename TFractFVGeom::SCVF& scvf = m_pFractGeo->scvf (ip);
  
  //  Compute K (rho * g) / mu
    VecScale (DarcyVel, m_Gravity, m_density * m_imFractPermeability[0] / m_viscosity);
 
  //  Compute flux
    number flux = VecDot (DarcyVel, scvf.normal ()) * m_density * half_fr_width;
 
  //  Add contribution to flow equation
    b (_P_, m_innerSideCo [scvf.from ()]) -= flux;
    b (_P_, m_innerSideCo [scvf.to ()])   += flux;
  }
  
//  mass sources
}
 
//  register assembling functions
 
template<typename TDomain>
template<typename TElem, typename TBulkFVGeom> // the element to register for and its full-dim. FV geometry
void FractPF_FV1<TDomain>::register_loc_discr_func_with_fvgeom ()
{
  static const ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
  
  this->clear_add_fct(id);
  
  this->set_prep_elem_loop_fct(id, & this_type::template prepare_element_loop<TElem, TBulkFVGeom>);
  this->set_prep_elem_fct   (id, & this_type::template prepare_element<TElem, TBulkFVGeom>);
  this->set_fsh_elem_loop_fct (id, & this_type::template finish_element_loop<TElem, TBulkFVGeom>);
  this->set_add_jac_A_elem_fct(id, & this_type::template ass_JA_elem<TElem, TBulkFVGeom>);
  this->set_add_jac_M_elem_fct(id, & this_type::template ass_JM_elem<TElem, TBulkFVGeom>);
  this->set_add_def_A_elem_fct(id, & this_type::template ass_dA_elem<TElem, TBulkFVGeom>);
  this->set_add_def_M_elem_fct(id, & this_type::template ass_dM_elem<TElem, TBulkFVGeom>);
  this->set_add_rhs_elem_fct  (id, & this_type::template ass_rhs_elem<TElem, TBulkFVGeom>);
}
 
template<typename TDomain>
template<typename TElem> // the element to register for
void FractPF_FV1<TDomain>::register_loc_discr_func ()
{
  if (! m_bCondensedFV)
    this->register_loc_discr_func_with_fvgeom<TElem, FV1Geometry<TElem, dim> > ();
  else
    this->register_loc_discr_func_with_fvgeom<TElem, FV1CondensedGeometry<TElem, dim> > ();
}
 
//  register assembling functions & parameters
 
template<typename TDomain>
void FractPF_FV1<TDomain>::register_discr_and_param (bool condensed)
{
//  initilize the 'condensed scvf ip' flag
  m_bCondensedFV = condensed;
  
//  register imports
  
  this->register_import (m_imPermeability);
  
  this->register_import (m_imAperture);
  this->register_import (m_imFractPermeability);
  
//  register the local assembling functions
  boost::mpl::for_each<AssembleElemList> (RegisterLocalDiscr (this));
}
 
template <typename TDomain>
FractPF_FV1<TDomain>::FractPF_FV1
(
  const char* functions, 
  const char* subsets, 
  bool condensed_FV 
)
: IElemDisc<TDomain> (functions, subsets),
  m_imPermeability (false),
  m_imAperture (false), m_imFractPermeability (false)
{
//  check number of functions
  if (this->num_fct () != 1)
    UG_THROW ("Wrong number of functions: The ElemDisc 'FractPF_FV1'"
          " needs exactly 1 symbolic function (for the pressure).");
 
//  register functions and data
  this->register_discr_and_param (condensed_FV);
};
 
template <typename TDomain>
FractPF_FV1<TDomain>::FractPF_FV1
(
  const char* functions, 
  const char* subsets 
)
: IElemDisc<TDomain> (functions, subsets),
  m_imPermeability (false),
  m_imAperture (false), m_imFractPermeability (false)
{
//  check number of functions
  if (this->num_fct () != 1)
    UG_THROW ("Wrong number of functions: The ElemDisc 'FractPF_FV1'"
          " needs exactly 1 symbolic function (for the pressure).");
 
//  register functions and data
  this->register_discr_and_param (false);
};
 
} // namespace d3f
} // end namespace ug
 
/* End of File */