inout_bnd_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
 */
 
 
namespace ug {
namespace d3f {
  
/*----- Class OutFlowBC -----*/
  
template <typename TDomain, typename TAlgebra>
void OutFlowBC<TDomain, TAlgebra>::add
(
  std::string ss_names, 
  std::string fct_names 
)
{
  std::vector<std::string> v_ss_names;
  std::vector<std::string> v_fct_names;
  
  TokenizeTrimString (ss_names, v_ss_names);
  for (size_t i = 0; i < v_ss_names.size(); i++)
  {
    v_fct_names.clear();
    TokenizeTrimString (fct_names, v_fct_names);
    if (v_fct_names.empty ()) continue;
    std::vector<std::string> & face_fct = m_mStoF[v_ss_names[i]];
    for (size_t j = 0; j < v_fct_names.size(); j++)
      face_fct.push_back(v_fct_names[j]);
  }
}
 
template <typename TDomain, typename TAlgebra>
void OutFlowBC<TDomain, TAlgebra>::scale
(
  std::string fct_names, 
  number scaling 
)
{
  std::vector<std::string> v_fct_names;
  TokenizeTrimString (fct_names, v_fct_names);
  for (size_t i = 0; i < v_fct_names.size (); i++)
    m_mFtoM [v_fct_names [i]] = scaling;
}
  
template <typename TDomain, typename TAlgebra>
void OutFlowBC<TDomain, TAlgebra>::names_to_ids_in_scaling
(
  const FunctionPattern * fct_pattern, 
  std::vector<number> & v_scaling 
) const
{
  typedef fct_name_scaling_map_type::const_iterator scaling_iter_type;
  
//  set the default scaling
  v_scaling.resize (fct_pattern->num_fct ());
  for (size_t i = 0; i < v_scaling.size (); i++)
    v_scaling[i] = 1;
  
//  loop the specifications and extract the scaling
  for (scaling_iter_type iter = m_mFtoM.begin (); iter != m_mFtoM.end (); ++iter)
  {
    const char * fct_name = iter->first.c_str ();
    number fct_scale = iter->second;
    
    for (size_t i = 0; i < v_scaling.size (); i++)
      if (strcmp (fct_name, fct_pattern->name (i)) == 0)
        v_scaling [i] = fct_scale;
  }
}
 
template <typename TDomain, typename TAlgebra>
void OutFlowBC<TDomain, TAlgebra>::adjust_jacobian
(
  matrix_type & J,
  const vector_type & u,
  ConstSmartPtr<DoFDistribution> dd,
  int type,
  number time,
  ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
  const number s_a0
)
{
  if (m_spVolStabData.invalid())
    UG_THROW ("OutFlowBC: No volume stabilization specified!");
  
  typedef bc_data_map_type::iterator t_face_iter;
  
  std::vector<DoFIndex> multInd;
  SubsetGroup ssGrp (base_type::m_spApproxSpace->subset_handler ());
  FunctionGroup fctGrp (base_type::m_spApproxSpace->function_pattern ());
  
//  Get the scaling of the equations
  std::vector<number> vScaling;
  names_to_ids_in_scaling (base_type::m_spApproxSpace->function_pattern().get (), vScaling);
 
//  Loop the entries of the map
  for (t_face_iter face = m_mStoF.begin(); face != m_mStoF.end(); ++face)
  {
    ssGrp.clear ();
    ssGrp.add (face->first);
    
    fctGrp.clear ();
    fctGrp.add (face->second);
    
  //  Loop over the subsets
    for (size_t i = 0; i < ssGrp.size (); i++)
    {
      int si = ssGrp [i];
      
    //  Loop the edges
      t_vert_iterator iterEnd = dd->end<Vertex> (si);
      for (t_vert_iterator iter = dd->begin<Vertex> (si); iter != iterEnd; iter++)
      {
        Vertex * vrt = *iter;
        MathVector<2> & flux_parts = m_spVolStabData->value(vrt);
        number flux = s_a0 * flux_parts[0] + flux_parts[1];
        
      //  Loop functions
        for(size_t j = 0; j < fctGrp.size (); j++)
        {
          size_t fct = fctGrp [j];
          
        //  Get the multiindices
          if (dd->inner_dof_indices (vrt, fct, multInd) != 1)
            UG_THROW ("OutFlowBC: More than one DoF per vertex. Not the Lagrange element?");
    
        //  Add the flux to the diagonal
          DoFRef (J, multInd[0], multInd[0]) += flux * vScaling[fct];
        }
      }
    }
  }
}
 
template <typename TDomain, typename TAlgebra>
void OutFlowBC<TDomain, TAlgebra>::adjust_defect
(
  vector_type& d,
  const vector_type& u,
  ConstSmartPtr<DoFDistribution> dd,
  int type,
  number time,
  ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
  const std::vector<number> * vScaleMass,
  const std::vector<number> * vScaleStiff
)
{
  if (vScaleMass == NULL || vScaleStiff == NULL)
    UG_THROW ("OutFlowBC: This BC can be used for instationary problems only!");
  if (m_spVolStabData.invalid())
    UG_THROW ("OutFlowBC: No volume stabilization specified!");
  
  typedef bc_data_map_type::iterator t_face_iter;
  
  std::vector<DoFIndex> multInd;
  SubsetGroup ssGrp (base_type::m_spApproxSpace->subset_handler ());
  FunctionGroup fctGrp (base_type::m_spApproxSpace->function_pattern ());
  
//  Get the scaling of the equations
  std::vector<number> vScaling;
  names_to_ids_in_scaling (base_type::m_spApproxSpace->function_pattern().get (), vScaling);
 
//  Loop the entries of the map
  for (t_face_iter face = m_mStoF.begin(); face != m_mStoF.end(); ++face)
  {
    ssGrp.clear ();
    ssGrp.add (face->first);
    
    fctGrp.clear ();
    fctGrp.add (face->second);
    
  //  Loop over the subsets
    for (size_t i = 0; i < ssGrp.size (); i++)
    {
      int si = ssGrp [i];
      
    //  Loop the edges
      t_vert_iterator iterEnd = dd->end<Vertex> (si);
      for (t_vert_iterator iter = dd->begin<Vertex> (si); iter != iterEnd; iter++)
      {
        Vertex * vrt = *iter;
        MathVector<2> & flux_part = m_spVolStabData->value(vrt);
        number flux = (* vScaleStiff)[0] * flux_part[0] + (* vScaleMass)[0] * flux_part[1];
        
      //  Loop functions
        for(size_t j = 0; j < fctGrp.size (); j++)
        {
          size_t fct = fctGrp [j];
          
        //  Get the multiindices
          if (dd->inner_dof_indices (vrt, fct, multInd) != 1)
            UG_THROW ("OutFlowBC: More than one DoF per vertex. Not the Lagrange element?");
    
        //  Take into account the outflow flux
          DoFRef (d, multInd[0]) += flux * vScaling[fct] * DoFRef (u, multInd[0]);
        }
      }
    }
  }
}
 
/*----- Class InOutFlowBC -----*/
 
template <typename TDomain, typename TAlgebra>
void InOutFlowBC<TDomain, TAlgebra>::add_spec
(
  std::string ss_names, 
  std::string fct_names, 
  number val, 
  SmartPtr<UserData<number, dim> > func 
)
{
  std::vector<std::string> v_ss_names;
  std::vector<std::string> v_fct_names;
  
  TokenizeTrimString (ss_names, v_ss_names);
  for (size_t i = 0; i < v_ss_names.size(); i++)
  {
    v_fct_names.clear();
    TokenizeTrimString (fct_names, v_fct_names);
    if (v_fct_names.empty ()) continue;
    std::map<std::string, t_inDirichlet> & face_fct = m_mStoF[v_ss_names[i]];
    for (size_t j = 0; j < v_fct_names.size(); j++)
    {
      t_inDirichlet & data = face_fct[v_fct_names[j]];
      data.val = val;
      data.func = func;
    }
  }
}
 
template <typename TDomain, typename TAlgebra>
void InOutFlowBC<TDomain, TAlgebra>::scale
(
  std::string fct_names, 
  number scaling 
)
{
  std::vector<std::string> v_fct_names;
  TokenizeTrimString (fct_names, v_fct_names);
  for (size_t i = 0; i < v_fct_names.size (); i++)
    m_mFtoM [v_fct_names [i]] = scaling;
}
  
template <typename TDomain, typename TAlgebra>
void InOutFlowBC<TDomain, TAlgebra>::names_to_ids_in_scaling
(
  const FunctionPattern * fct_pattern, 
  std::vector<number> & v_scaling 
) const
{
  typedef fct_name_scaling_map_type::const_iterator scaling_iter_type;
  
//  set the default scaling
  v_scaling.resize (fct_pattern->num_fct ());
  for (size_t i = 0; i < v_scaling.size (); i++)
    v_scaling[i] = 1;
  
//  loop the specifications and extract the scaling
  for (scaling_iter_type iter = m_mFtoM.begin (); iter != m_mFtoM.end (); ++iter)
  {
    const char * fct_name = iter->first.c_str ();
    number fct_scale = iter->second;
    
    for (size_t i = 0; i < v_scaling.size (); i++)
      if (strcmp (fct_name, fct_pattern->name (i)) == 0)
        v_scaling [i] = fct_scale;
  }
}
 
template <typename TDomain, typename TAlgebra>
void InOutFlowBC<TDomain, TAlgebra>::adjust_jacobian
(
  matrix_type & J,
  const vector_type & u,
  ConstSmartPtr<DoFDistribution> dd,
  int type,
  number time,
  ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
  const number s_a0
)
{
  if (m_spVolStabData.invalid())
    UG_THROW ("InOutFlowBC: No volume stabilization specified!");
  
  typedef typename bc_data_map_type::iterator t_face_iter;
  typedef typename std::map<std::string, t_inDirichlet>::iterator t_func_iter;
  
  std::vector<DoFIndex> multInd;
  SubsetGroup ssGrp (base_type::m_spApproxSpace->subset_handler ());
  FunctionGroup fctGrp (base_type::m_spApproxSpace->function_pattern ());
  
//  Get the scaling of the equations
  std::vector<number> vScaling;
  names_to_ids_in_scaling (base_type::m_spApproxSpace->function_pattern().get (), vScaling);
 
//  Loop the entries of the map
  for (t_face_iter face = m_mStoF.begin(); face != m_mStoF.end(); ++face)
  {
    ssGrp.clear ();
    ssGrp.add (face->first);
    
    fctGrp.clear ();
    t_func_iter func_end = face->second.end();
    for (t_func_iter func_iter = face->second.begin(); func_iter != func_end; ++func_iter)
      fctGrp.add (func_iter->first);
    
  //  Loop over the subsets
    for (size_t i = 0; i < ssGrp.size (); i++)
    {
      int si = ssGrp [i];
      
    //  Loop the edges
      t_vert_iterator iterEnd = dd->end<Vertex> (si);
      for (t_vert_iterator iter = dd->begin<Vertex> (si); iter != iterEnd; iter++)
      {
        Vertex * vrt = *iter;
        MathVector<2> & flux_parts = m_spVolStabData->value(vrt);
        number flux = s_a0 * flux_parts[0] + flux_parts[1];
        
      //  Loop functions
        for(size_t j = 0; j < fctGrp.size (); j++)
        {
          size_t fct = fctGrp [j];
          
        //  Get the multiindices
          if (dd->inner_dof_indices (vrt, fct, multInd) != 1)
            UG_THROW ("InOutFlowBC: More than one DoF per vertex. Not the Lagrange element?");
        
          if (flux > 0)
          {
          //  Add the flux to the diagonal
            DoFRef (J, multInd[0], multInd[0]) += flux * vScaling[fct];
          //  Mark the vertex as an outflow vertex
            m_spVolStabData->set_flag (vrt, true);
          }
          else
          {
          //  Dirichlet: Set the identity matrix
            SetDirichletRow (J, multInd[0]);
          //  Mark the vertex as a Dirichlet vertex
            m_spVolStabData->set_flag (vrt, false);
          }
        }
      }
    }
  }
}
 
template <typename TDomain, typename TAlgebra>
void InOutFlowBC<TDomain, TAlgebra>::adjust_defect
(
  vector_type& d,
  const vector_type& u,
  ConstSmartPtr<DoFDistribution> dd,
  int type,
  number time,
  ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
  const std::vector<number> * vScaleMass,
  const std::vector<number> * vScaleStiff
)
{
  if (vScaleMass == NULL || vScaleStiff == NULL)
    UG_THROW ("InOutFlowBC: This BC can be used for instationary problems only!");
  if (m_spVolStabData.invalid())
    UG_THROW ("InOutFlowBC: No volume stabilization specified!");
  
  typedef typename bc_data_map_type::iterator t_face_iter;
  typedef typename std::map<std::string, t_inDirichlet>::iterator t_func_iter;
  
  std::vector<DoFIndex> multInd;
  SubsetGroup ssGrp (base_type::m_spApproxSpace->subset_handler ());
  FunctionGroup fctGrp (base_type::m_spApproxSpace->function_pattern ());
  
//  Get the scaling of the equations
  std::vector<number> vScaling;
  names_to_ids_in_scaling (base_type::m_spApproxSpace->function_pattern().get (), vScaling);
 
//  Loop the entries of the map
  for (t_face_iter face = m_mStoF.begin(); face != m_mStoF.end(); ++face)
  {
    ssGrp.clear ();
    ssGrp.add (face->first);
    
    fctGrp.clear ();
    t_func_iter func_end = face->second.end();
    for (t_func_iter func_iter = face->second.begin(); func_iter != func_end; ++func_iter)
      fctGrp.add (func_iter->first);
    
  //  Loop over the subsets
    for (size_t i = 0; i < ssGrp.size (); i++)
    {
      int si = ssGrp [i];
      
    //  Loop the edges
      t_vert_iterator iterEnd = dd->end<Vertex> (si);
      for (t_vert_iterator iter = dd->begin<Vertex> (si); iter != iterEnd; iter++)
      {
        Vertex * vrt = *iter;
        MathVector<2> & flux_part = m_spVolStabData->value(vrt);
        number flux = (* vScaleStiff)[0] * flux_part[0] + (* vScaleMass)[0] * flux_part[1];
        
      //  Loop functions
        for(size_t j = 0; j < fctGrp.size (); j++)
        {
          size_t fct = fctGrp [j];
          
        //  Get the multiindices
          if (dd->inner_dof_indices (vrt, fct, multInd) != 1)
            UG_THROW ("InOutFlowBC: More than one DoF per vertex. Not the Lagrange element?");
          
          if (flux > 0)
          {
          //  Take into account the outflow flux
            DoFRef (d, multInd[0]) += flux * vScaling[fct] * DoFRef (u, multInd[0]);
          //  Mark the vertex as an outflow vertex
            m_spVolStabData->set_flag (vrt, true);
          }
          else
          {
          //  Dirichlet
            DoFRef (d, multInd[0]) = 0;
          //  Mark the vertex as a Dirichlet vertex
            m_spVolStabData->set_flag (vrt, false);
          }
        }
      }
    }
  }
}
 
template <typename TDomain, typename TAlgebra>
void InOutFlowBC<TDomain, TAlgebra>::adjust_solution
(
  vector_type & u,
  ConstSmartPtr<DoFDistribution> dd,
  int type,
  number time
)
{
  if (m_spVolStabData.invalid())
    UG_THROW ("InOutFlowBC: No volume stabilization specified!");
  
  typedef typename bc_data_map_type::iterator t_face_iter;
  typedef typename std::map<std::string, t_inDirichlet>::iterator t_func_iter;
  
  std::vector<DoFIndex> multInd;
  SubsetGroup ssGrp (base_type::m_spApproxSpace->subset_handler ());
  ConstSmartPtr<FunctionPattern> fctPattern = base_type::m_spApproxSpace->function_pattern ();
  typename domain_type::position_accessor_type aaPos = base_type::m_spApproxSpace->domain()->position_accessor();
  
//  Loop the entries of the map
  for (t_face_iter face = m_mStoF.begin(); face != m_mStoF.end(); ++face)
  {
    ssGrp.clear ();
    ssGrp.add (face->first);
    
  //  Loop over the subsets
    for (size_t i = 0; i < ssGrp.size (); i++)
    {
      int si = ssGrp [i];
      
    //  Loop the edges
      t_vert_iterator iterEnd = dd->end<Vertex> (si);
      for (t_vert_iterator iter = dd->begin<Vertex> (si); iter != iterEnd; iter++)
      {
        Vertex * vrt = *iter;
        if (m_spVolStabData->flag (vrt)) continue; // this is outflow, we do nothing here
        
      //  Loop functions
        t_func_iter func_end = face->second.end();
        for (t_func_iter func_iter = face->second.begin(); func_iter != func_end; ++func_iter)
        {
        //  Get the index of the function
          size_t fct_id = fctPattern->fct_id_by_name (func_iter->first.c_str ());
          
        //  Get the multiindices
          if (dd->inner_dof_indices (vrt, fct_id, multInd) != 1)
            UG_THROW ("InOutFlowBC: More than one DoF per vertex. Not the Lagrange element?");
          
        //  Get the Dirichlet value for this vertex
          t_inDirichlet & inflow_data = func_iter->second;
          if (inflow_data.func.valid ())
            (* inflow_data.func) (inflow_data.val, aaPos[vrt], time, si);
        
        //  Set the value
          DoFRef (u, multInd[0]) = inflow_data.val;
          /*TODO: These values must depend on the function! */
        }
      }
    }
  }
}
 
} // namespace d3f
} // end namespace ug
 
/* End of File */