diff --git a/CMakeLists.txt b/CMakeLists.txt index 2501d6c..89e759c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -38,6 +38,8 @@ set(pluginName ConvectionDiffusion) set(SOURCES convection_diffusion_base.cpp fv1/convection_diffusion_fv1.cpp + bnd/outflow_base.cpp + fv1/bnd/outflow_fv1.cpp fe/convection_diffusion_fe.cpp fe/convection_diffusion_stab_fe.cpp fvcr/convection_diffusion_fvcr.cpp @@ -98,4 +100,4 @@ if(USE_PYBIND11) # CPack set_target_properties(pyconvectiondiffusion PROPERTIES INSTALL_RPATH "$ORIGIN/..:$ORIGIN/../../../lib") install(TARGETS pyconvectiondiffusion LIBRARY DESTINATION ug4py COMPONENT pymodules) -endif(USE_PYBIND11) \ No newline at end of file +endif(USE_PYBIND11) diff --git a/bnd/outflow_base.cpp b/bnd/outflow_base.cpp new file mode 100644 index 0000000..51fc734 --- /dev/null +++ b/bnd/outflow_base.cpp @@ -0,0 +1,136 @@ +/* + * Copyright (c) 2012-2013: G-CSC, Goethe University Frankfurt + * Authors: Daniel Gonzalez + * Based on the modules by Dmitry Logashenko, Andreas Vogel + * + * This file is part of UG4. + * + * UG4 is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License version 3 (as published by the + * Free Software Foundation) with the following additional attribution + * requirements (according to LGPL/GPL v3 §7): + * + * (1) The following notice must be displayed in the Appropriate Legal Notices + * of covered and combined works: "Based on UG4 (www.ug4.org/license)". + * + * (2) The following notice must be displayed at a prominent place in the + * terminal output of covered works: "Based on UG4 (www.ug4.org/license)". + * + * (3) The following bibliography is recommended for citation and must be + * preserved in all covered files: + * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively + * parallel geometric multigrid solver on hierarchically distributed grids. + * Computing and visualization in science 16, 4 (2013), 151-164" + * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel + * flexible software system for simulating pde based models on high performance + * computers. Computing and visualization in science 16, 4 (2013), 165-179" + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + */ + + +#include "outflow_base.h" + +#include "common/util/provider.h" +#include "lib_disc/spatial_disc/disc_util/fv1_geom.h" + +namespace ug{ +namespace ConvectionDiffusionPlugin{ + +/** + * converts the subset names where the BC is imposed to the corresponding subset + * indices (i.e. m_vScheduledBndSubSets -> m_vBndSubSetIndex): + */ +template +void ConvectionDiffusionOutflowBase::extract_scheduled_data() +{ +// clear all extracted data + m_vBndSubSetIndex.clear(); + +// loop all scheduled subsets + for(size_t i = 0; i < m_vScheduledBndSubSets.size(); ++i) + { + // create Subset Group + SubsetGroup subsetGroup; + + // convert strings + try{ + subsetGroup = this->approx_space()->subset_grp_by_name(m_vScheduledBndSubSets[i].c_str()); + }UG_CATCH_THROW("'OutflowBase:extract_scheduled_data':" + " Subsets '" <function_pattern()->subset_handler(); + + // loop subsets + for(size_t si = 0; si < subsetGroup.size(); ++si) + { + // get subset index + const int subsetIndex = subsetGroup[si]; + + // check that subsetIndex is valid + if(subsetIndex < 0 || subsetIndex >= rSH.num_subsets()) + { + UG_LOG("ERROR in 'OutflowBase:extract_scheduled_data':" + " Invalid subset Index " << subsetIndex << + ". (Valid is 0, .. , " << rSH.num_subsets() <<").\n"); + return; + } + + // save the index + m_vBndSubSetIndex.push_back(subsetIndex); + } + } +} + +/** + * The add method for the boundary subsets: + */ +template +void ConvectionDiffusionOutflowBase::add +( + const char* subsets // string with the ','-separated names of the subsets +) +{ + m_vScheduledBndSubSets.push_back(subsets); +} + +//////////////////////////////////////////////////////////////////////////////// +// Constructor +//////////////////////////////////////////////////////////////////////////////// + +template +ConvectionDiffusionOutflowBase:: +ConvectionDiffusionOutflowBase(SmartPtr< ConvectionDiffusionBase > spMaster) + : IElemDisc(spMaster->symb_fcts(), spMaster->symb_subsets()), + m_spMaster (spMaster) +{ +// check number of functions + if(this->num_fct() != 1) + UG_THROW("Wrong number of functions: The ElemDisc 'ConvectionDiffusion'" + " needs exactly "<<1<<" symbolic function."); + +// yet no boundary subsets + m_vBndSubSetIndex.clear (); +} + +//////////////////////////////////////////////////////////////////////////////// +// explicit template instantiations +//////////////////////////////////////////////////////////////////////////////// + +#ifdef UG_DIM_1 +template class ConvectionDiffusionOutflowBase; +#endif +#ifdef UG_DIM_2 +template class ConvectionDiffusionOutflowBase; +#endif +#ifdef UG_DIM_3 +template class ConvectionDiffusionOutflowBase; +#endif + +} // namespace ConvectionDiffusionPlugin +} // namespace ug diff --git a/bnd/outflow_base.h b/bnd/outflow_base.h new file mode 100644 index 0000000..c83439a --- /dev/null +++ b/bnd/outflow_base.h @@ -0,0 +1,115 @@ +/* + * Copyright (c) 2012-2013: G-CSC, Goethe University Frankfurt + * Authors: Daniel Gonzalez + * Based on the modules by Dmitry Logashenko, Andreas Vogel + * + * This file is part of UG4. + * + * UG4 is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License version 3 (as published by the + * Free Software Foundation) with the following additional attribution + * requirements (according to LGPL/GPL v3 §7): + * + * (1) The following notice must be displayed in the Appropriate Legal Notices + * of covered and combined works: "Based on UG4 (www.ug4.org/license)". + * + * (2) The following notice must be displayed at a prominent place in the + * terminal output of covered works: "Based on UG4 (www.ug4.org/license)". + * + * (3) The following bibliography is recommended for citation and must be + * preserved in all covered files: + * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively + * parallel geometric multigrid solver on hierarchically distributed grids. + * Computing and visualization in science 16, 4 (2013), 151-164" + * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel + * flexible software system for simulating pde based models on high performance + * computers. Computing and visualization in science 16, 4 (2013), 165-179" + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + */ + +#ifndef __H__UG__PLUGINS__CONVECTION_DIFFUSION__CONVECTION_DIFFUSION_FV1__BND__OUTFLOW_BASE__ +#define __H__UG__PLUGINS__CONVECTION_DIFFUSION__CONVECTION_DIFFUSION_FV1__BND__OUTFLOW_BASE__ +// other ug4 modules +#include "common/common.h" +#include "lib_grid/lg_base.h" + +// library intern headers +#include "lib_disc/spatial_disc/elem_disc/elem_disc_interface.h" +#include "lib_disc/spatial_disc/user_data/data_export.h" +#include "lib_disc/spatial_disc/user_data/data_import.h" + +#include "../convection_diffusion_base.h" + +namespace ug{ +namespace ConvectionDiffusionPlugin{ + + +/// \ingroup lib_disc_elem_disc +/// @{ + +/// The zero-stress (neutral) outflow boundary condition for the incompressible NS equation +/** + * This class implements the so-called neutral boundary condition for outflow + * boundaries. Note that this class can be used only with the stabilized + * vertex-centered discretization of the Navier-Stokes equations. + * + * This boundary condition imposes two equations on the unknown functions + * at the outflow boundary \f$ F \f$: + *
    + *
  • \f$ \int_F p ds = 0 \f$, and + *
  • \f$ \sigma \vec{n} \cdot \vec{n} = 0 \f$ on \f$ F \f$. + *
+ * where \f$ \vec{n} \f$ is the outer normal at \f$ F \f$ and + * \f$ \sigma = \mu (\nabla \vec{u} + (\nabla \vec{u})^T) \f$ the stress tensor. + */ +template< typename TDomain> +class ConvectionDiffusionOutflowBase + : public IElemDisc +{ + protected: + /// Base class type + typedef IElemDisc base_type; + + /// own type + typedef ConvectionDiffusionOutflowBase this_type; + + public: + /// World dimension + static const int dim = base_type::dim; + + public: + /// Constructor (setting default values) + ConvectionDiffusionOutflowBase(SmartPtr< ConvectionDiffusionBase > spMaster); + + /// adds a boundary segment + void add(const char* subsets); + + protected: + /// sets the kinematic viscosity + virtual void set_velocity(SmartPtr, dim> > data) = 0; + + public: + /// returns if local time series is needed + virtual bool requests_local_time_series() {return true;} + + protected: + /// The master discretization: + SmartPtr< ConvectionDiffusionBase > m_spMaster; + + /// The boundary subsets: + std::vector m_vScheduledBndSubSets; // names + std::vector m_vBndSubSetIndex; // indices + + void extract_scheduled_data(); // convert m_vScheduledBndSubSets -> m_vBndSubSetIndex +}; + +/// @} + +} // namespace NavierStokes +} // end namespace ug + +#endif /*_H__UG__PLUGINS__CONVECTION_DIFFUSION__CONVECTION_DIFFUSION_FV1__BND__OUTFLOW_BASE__*/ diff --git a/convection_diffusion_base.h b/convection_diffusion_base.h index 535beef..3725bd0 100644 --- a/convection_diffusion_base.h +++ b/convection_diffusion_base.h @@ -118,6 +118,8 @@ class ConvectionDiffusionBase void set_velocity(LuaFunctionHandle fct); #endif /// \} + /// /// returns velocity + virtual SmartPtr, dim> > velocity() = 0; /// sets the flux /** diff --git a/convection_diffusion_plugin.cpp b/convection_diffusion_plugin.cpp index a1cba0b..fc39090 100644 --- a/convection_diffusion_plugin.cpp +++ b/convection_diffusion_plugin.cpp @@ -39,6 +39,7 @@ #include "convection_diffusion_base.h" #include "convection_diffusion_sss.h" #include "fv1/convection_diffusion_fv1.h" +#include "fv1/bnd/outflow_fv1.h" #include "fe/convection_diffusion_fe.h" #include "fe/convection_diffusion_stab_fe.h" #include "fvcr/convection_diffusion_fvcr.h" @@ -193,6 +194,26 @@ static void Domain(TRegistry& reg, string grp) .set_construct_as_smart_pointer(true); reg.add_class_to_group(name, "ConvectionDiffusionFV1", tag); } + + // Convection Diffusion FV1 Outflow boundary condition base + { + typedef ConvectionDiffusionOutflowBase T; + typedef IElemDisc TBase; + string name = string("ConvectionDiffusionOutflowBase").append(suffix); + reg.template add_class_(name, grp) + .add_method("add", &T::add, "", "Subset(s)"); + reg.add_class_to_group(name, "ConvectionDiffusionOutflowBase", tag); + } +// Convection Diffusion FV1 Outflow boundary condition + { + typedef ConvectionDiffusionOutflowFV1 T; + typedef ConvectionDiffusionOutflowBase TBase; + string name = string("ConvectionDiffusionOutflowFV1").append(suffix); + reg.template add_class_(name, grp) + .template add_constructor >)>("MasterDisc") + .set_construct_as_smart_pointer(true); + reg.add_class_to_group(name, "ConvectionDiffusionOutflowFV1", tag); + } // Convection Diffusion FE { diff --git a/fe/convection_diffusion_fe.h b/fe/convection_diffusion_fe.h index bd4d4c9..d3672b3 100644 --- a/fe/convection_diffusion_fe.h +++ b/fe/convection_diffusion_fe.h @@ -76,7 +76,9 @@ class ConvectionDiffusionFE : public ConvectionDiffusionBase /// sets the quad order void set_quad_order(size_t order); - + + /// returns velocity + SmartPtr, dim> > velocity() {return m_imVelocity.user_data ();} private: /// prepares the loop over all elements /** diff --git a/fractfv1/convection_diffusion_fractfv1.h b/fractfv1/convection_diffusion_fractfv1.h index 3175f36..70aa522 100644 --- a/fractfv1/convection_diffusion_fractfv1.h +++ b/fractfv1/convection_diffusion_fractfv1.h @@ -140,6 +140,9 @@ class ConvectionDiffusionFractFV1 : public ConvectionDiffusionBase } #endif + /// returns velocity + SmartPtr, dim> > velocity() {return m_imVelocity.user_data ();} + void set_ortho_velocity(SmartPtr > user) { m_imOrthoVelocity.set_data(user); diff --git a/fv/convection_diffusion_fv.h b/fv/convection_diffusion_fv.h index e8defcc..c2d21fa 100644 --- a/fv/convection_diffusion_fv.h +++ b/fv/convection_diffusion_fv.h @@ -225,6 +225,8 @@ class ConvectionDiffusionFV : public ConvectionDiffusionBase /// sets the quad order void set_quad_order(size_t order); + /// returns velocity + SmartPtr, dim> > velocity() {return m_imVelocity.user_data ();} protected: /// current shape function set diff --git a/fv1/bnd/outflow_fv1.cpp b/fv1/bnd/outflow_fv1.cpp new file mode 100644 index 0000000..a8b12e1 --- /dev/null +++ b/fv1/bnd/outflow_fv1.cpp @@ -0,0 +1,422 @@ +/* + * Copyright (c) 2012-2014: G-CSC, Goethe University Frankfurt + * Authors: Daniel Gonzalez + * Based on the modules by Dmitry Logashenko, Andreas Vogel + * + * This file is part of UG4. + * + * UG4 is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License version 3 (as published by the + * Free Software Foundation) with the following additional attribution + * requirements (according to LGPL/GPL v3 §7): + * + * (1) The following notice must be displayed in the Appropriate Legal Notices + * of covered and combined works: "Based on UG4 (www.ug4.org/license)". + * + * (2) The following notice must be displayed at a prominent place in the + * terminal output of covered works: "Based on UG4 (www.ug4.org/license)". + * + * (3) The following bibliography is recommended for citation and must be + * preserved in all covered files: + * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively + * parallel geometric multigrid solver on hierarchically distributed grids. + * Computing and visualization in science 16, 4 (2013), 151-164" + * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel + * flexible software system for simulating pde based models on high performance + * computers. Computing and visualization in science 16, 4 (2013), 165-179" + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + */ + + +#include "outflow_fv1.h" + +#include "lib_disc/spatial_disc/disc_util/fv1_geom.h" +#include "lib_disc/spatial_disc/disc_util/geom_provider.h" + +namespace ug{ +namespace ConvectionDiffusionPlugin{ + +//////////////////////////////////////////////////////////////////////////////// +// Constructor - set default values +//////////////////////////////////////////////////////////////////////////////// + +template +ConvectionDiffusionOutflowFV1:: +ConvectionDiffusionOutflowFV1(SmartPtr< ConvectionDiffusionBase > spMaster) +: ConvectionDiffusionOutflowBase(spMaster) +{ + m_imVelocity.set_comp_lin_defect(false); + +// register imports + this->register_import(m_imVelocity); + + +// initialize the imports from the master discretization + set_velocity(spMaster->velocity ()); + + // update assemble functions + register_all_funcs(false); +}; + + +template +void ConvectionDiffusionOutflowFV1:: +prepare_setting(const std::vector& vLfeID, bool bNonRegularGrid) +{ + if(bNonRegularGrid) + UG_THROW("ConvectionDiffusion: only regular grid implemented."); + +// check number + if(vLfeID.size() != 1) + UG_THROW("ConvectionDiffusion: Need exactly "<<1<<" functions"); + + if(vLfeID[0].type() != LFEID::LAGRANGE || vLfeID[0].order() != 1) + UG_THROW("ConvectionDiffusion: 'fv1' expects Lagrange P1 trial space " + "for concentration."); + + + // update assemble functions + register_all_funcs(false); +} + +//////////////////////////////////////////////////////////////////////////////// +// assembling functions +//////////////////////////////////////////////////////////////////////////////// + +/** + * Prepares the element loop for a given element type: computes the FV-geo, ... + * Note that there are separate loops for every type of the grid elements. + */ +template +template +void ConvectionDiffusionOutflowFV1:: +prep_elem_loop(const ReferenceObjectID roid, const int si) +{ +// register subsetIndex at Geometry + static TFVGeom& geo = GeomProvider::get(); + +// Only first order implementation + if(!(TFVGeom::order == 1)) + UG_THROW("Only first order implementation, but other Finite Volume" + " Geometry set."); + +// check if velocity has been set + if(!m_imVelocity.data_given()) + UG_THROW("OutflowFV1::prep_elem_loop:" + " velocity has not been set, but is required.\n"); + +// extract indices of boundary + this->extract_scheduled_data(); + +// request the subset indices as boundary subset. This will force the +// creation of boundary subsets when calling geo.update + typename std::vector::const_iterator subsetIter; + for(subsetIter = m_vBndSubSetIndex.begin(); + subsetIter != m_vBndSubSetIndex.end(); ++subsetIter) + geo.add_boundary_subset(*subsetIter); +} + +/** + * Finalizes the element loop for a given element type. + */ +template +template +void ConvectionDiffusionOutflowFV1:: +fsh_elem_loop() +{ + static TFVGeom& geo = GeomProvider::get(); + +// remove the bnd subsets + typename std::vector::const_iterator subsetIter; + for(subsetIter = m_vBndSubSetIndex.begin(); + subsetIter != m_vBndSubSetIndex.end(); ++subsetIter) + geo.remove_boundary_subset(*subsetIter); +} + + +/** + * General initializations of a given grid element for the assembling. + */ +template +template +void ConvectionDiffusionOutflowFV1:: +prep_elem(const LocalVector& u, GridObject* elem, const ReferenceObjectID roid, const MathVector vCornerCoords[]) +{ +// Update Geometry for this element + static TFVGeom& geo = GeomProvider::get(); + try{ + geo.update(elem, vCornerCoords, &(this->subset_handler())); + } + UG_CATCH_THROW("OutflowFV1::prep_elem:" + " Cannot update Finite Volume Geometry."); + +// find and set the local and the global positions of the IPs for imports + typedef typename TFVGeom::BF BF; + typename std::vector::const_iterator subsetIter; + + m_vLocIP.clear(); m_vGloIP.clear(); + for(subsetIter = m_vBndSubSetIndex.begin(); + subsetIter != m_vBndSubSetIndex.end(); ++subsetIter) + { + const int bndSubset = *subsetIter; + if(geo.num_bf(bndSubset) == 0) continue; + const std::vector& vBF = geo.bf(bndSubset); + for(size_t i = 0; i < vBF.size(); ++i) + { + m_vLocIP.push_back(vBF[i].local_ip()); + m_vGloIP.push_back(vBF[i].global_ip()); + } + } + // REMARK: The loop above determines the ordering of the integration points: + // The "outer ordering" corresponds to the ordering of the subsets in + // m_vBndSubSetIndex, and "inside" of this ordering, the ip's are ordered + // according to the order of the boundary faces in the FV geometry structure. + + m_imVelocity.set_local_ips(&m_vLocIP[0], m_vLocIP.size()); + m_imVelocity.set_global_ips(&m_vGloIP[0], m_vGloIP.size()); + +} + + +/// Assembling of the convective flux in the Jacobian of the trasport eq. +template +template +void ConvectionDiffusionOutflowFV1:: +convective_flux_Jac +( + const size_t ip, // index of the integration point (for the velocity) + const BF& bf, // boundary face to assemble + LocalMatrix& J, // local Jacobian to update + const LocalVector& u, // local solution + const number& stdValue // concentration at ip +) +{ +// The convection velocity according to the current approximation: + number volumetric_flux = VecDot (m_imVelocity[ip], bf.normal ()) * 1; + +// We assume that there should be no inflow through the outflow boundary: + if (volumetric_flux < 0) + volumetric_flux = 0; + +// Add flux to local Jacobian + for(size_t sh = 0; sh < bf.num_sh(); ++sh) + { + number t = volumetric_flux * bf.shape(sh); + J(_C_, bf.node_id(), _C_, sh) += t; + } +} + +/// Assembling of the convective flux in the defect of the trasport eq. +template +template +void ConvectionDiffusionOutflowFV1:: +convective_flux_defect +( + const size_t ip, // index of the integration point (for the velocity) + const BF& bf, // boundary face to assemble + LocalVector& d, // local defect to update + const LocalVector& u, // local solution + const number& stdValue // concetration at ip +) +{ +// The volumetric flux at boundary: + number volumetric_flux = VecDot (m_imVelocity[ip], bf.normal ()); + +// We assume that there should be no inflow through the outflow boundary: + if (volumetric_flux < 0) + volumetric_flux = 0; + +// Add the flux to the defect: + + d(_C_, bf.node_id()) += volumetric_flux * stdValue; +} + + + +template +template +void ConvectionDiffusionOutflowFV1:: +add_jac_A_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]) +{ +// Only first order implementation + UG_ASSERT((TFVGeom::order == 1), "Only first order implemented."); + +// get finite volume geometry + static const TFVGeom& geo = GeomProvider::get(); + typedef typename TFVGeom::BF BF; + +// loop registered boundary segments + typename std::vector::const_iterator subsetIter; + size_t ip = 0; + for(subsetIter = m_vBndSubSetIndex.begin(); + subsetIter != m_vBndSubSetIndex.end(); ++subsetIter) + { + // get subset index corresponding to boundary + const int bndSubset = *subsetIter; + + // get the list of the ip's: + if(geo.num_bf(bndSubset) == 0) continue; + const std::vector& vBF = geo.bf(bndSubset); + + // loop the boundary faces + typename std::vector::const_iterator bf; + for(bf = vBF.begin(); bf != vBF.end(); ++bf, ++ip) + { + + number stdValue=0; + convective_flux_Jac (ip, *bf, J, u, stdValue); + + + + } + } +} + +template +template +void ConvectionDiffusionOutflowFV1:: +add_def_A_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]) +{ +// Only first order implemented + UG_ASSERT((TFVGeom::order == 1), "Only first order implemented."); + +// get finite volume geometry + static const TFVGeom& geo = GeomProvider::get(); + typedef typename TFVGeom::BF BF; + +// loop registered boundary segments + typename std::vector::const_iterator subsetIter; + size_t ip = 0; + for(subsetIter = m_vBndSubSetIndex.begin(); + subsetIter != m_vBndSubSetIndex.end(); ++subsetIter) + { + // get subset index corresponding to boundary + const int bndSubset = *subsetIter; + + // get the list of the ip's: + if(geo.num_bf(bndSubset) == 0) continue; + const std::vector& vBF = geo.bf(bndSubset); + + // loop the boundary faces + typename std::vector::const_iterator bf; + for(bf = vBF.begin(); bf != vBF.end(); ++bf, ++ip) + { + + // A. Compute concentration at ip + number stdValue=0; + for(size_t sh = 0; sh < bf->num_sh(); ++sh) + { + stdValue += u(_C_, sh) * bf->shape(sh); + } + + convective_flux_defect (ip, *bf, d, u, stdValue); + + + } + } +} + + + + + +//////////////////////////////////////////////////////////////////////////////// +// register assemble functions +//////////////////////////////////////////////////////////////////////////////// + +#ifdef UG_DIM_1 +template<> +void ConvectionDiffusionOutflowFV1:: +register_all_funcs(bool bHang) +{ +// switch assemble functions + if(!bHang) + { + register_func >(); + } + else + { + UG_THROW("OutflowFV1: Hanging Nodes not implemented.") + } +} +#endif + +#ifdef UG_DIM_2 +template<> +void ConvectionDiffusionOutflowFV1:: +register_all_funcs(bool bHang) +{ +// switch assemble functions + if(!bHang) + { + register_func >(); + register_func >(); + } + else + { + UG_THROW("OutflowFV1: Hanging Nodes not implemented.") + } +} +#endif + +#ifdef UG_DIM_3 +template<> +void ConvectionDiffusionOutflowFV1:: +register_all_funcs(bool bHang) +{ +// switch assemble functions + if(!bHang) + { + register_func >(); + register_func >(); + register_func >(); + register_func >(); + } + else + { + UG_THROW("OutflowFV1: Hanging Nodes not implemented.") + } +} +#endif + +template +template +void +ConvectionDiffusionOutflowFV1:: +register_func() +{ + ReferenceObjectID id = geometry_traits::REFERENCE_OBJECT_ID; + typedef this_type T; + + this->clear_add_fct(id); + this->set_prep_elem_loop_fct( id, &T::template prep_elem_loop); + this->set_prep_elem_fct( id, &T::template prep_elem); + this->set_fsh_elem_loop_fct( id, &T::template fsh_elem_loop); + this->set_add_jac_A_elem_fct( id, &T::template add_jac_A_elem); + this->set_add_jac_M_elem_fct( id, &T::template add_jac_M_elem); + this->set_add_def_A_elem_fct( id, &T::template add_def_A_elem); + this->set_add_def_M_elem_fct( id, &T::template add_def_M_elem); + this->set_add_rhs_elem_fct( id, &T::template add_rhs_elem); +} + + +//////////////////////////////////////////////////////////////////////////////// +// explicit template instantiations +//////////////////////////////////////////////////////////////////////////////// + +#ifdef UG_DIM_1 +template class ConvectionDiffusionOutflowFV1; +#endif +#ifdef UG_DIM_2 +template class ConvectionDiffusionOutflowFV1; +#endif +#ifdef UG_DIM_3 +template class ConvectionDiffusionOutflowFV1; +#endif + +} // namespace ConvectionDiffusionPlugin +} // namespace ug diff --git a/fv1/bnd/outflow_fv1.h b/fv1/bnd/outflow_fv1.h new file mode 100644 index 0000000..500d16d --- /dev/null +++ b/fv1/bnd/outflow_fv1.h @@ -0,0 +1,183 @@ +/* + * Copyright (c) 2012-2014: G-CSC, Goethe University Frankfurt + * Authors: Daniel Gonzalez + * Based on the modules by Dmitry Logashenko, Andreas Vogel + * + * This file is part of UG4. + * + * UG4 is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License version 3 (as published by the + * Free Software Foundation) with the following additional attribution + * requirements (according to LGPL/GPL v3 §7): + * + * (1) The following notice must be displayed in the Appropriate Legal Notices + * of covered and combined works: "Based on UG4 (www.ug4.org/license)". + * + * (2) The following notice must be displayed at a prominent place in the + * terminal output of covered works: "Based on UG4 (www.ug4.org/license)". + * + * (3) The following bibliography is recommended for citation and must be + * preserved in all covered files: + * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively + * parallel geometric multigrid solver on hierarchically distributed grids. + * Computing and visualization in science 16, 4 (2013), 151-164" + * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel + * flexible software system for simulating pde based models on high performance + * computers. Computing and visualization in science 16, 4 (2013), 165-179" + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + */ + +#ifndef __H__UG__PLUGINS__CONVECTION_DIFFUSION__CONVECTION_DIFFUSION_FV1__BND__OUTFLOW_FV1__ +#define __H__UG__PLUGINS__CONVECTION_DIFFUSION__CONVECTION_DIFFUSION_FV1__BND__OUTFLOW_FV1__ + + +// other ug4 modules +#include "common/common.h" +#include "lib_grid/lg_base.h" + +// library intern headers +#include "lib_disc/spatial_disc/elem_disc/elem_disc_interface.h" +#include "lib_disc/spatial_disc/user_data/data_export.h" +#include "lib_disc/spatial_disc/user_data/data_import.h" + +#include "../../bnd/outflow_base.h" + +namespace ug{ +namespace ConvectionDiffusionPlugin{ + + +/// \ingroup lib_disc_elem_disc +/// @{ + +/// Outflow boundary condition for the Transport equation +/** + * This class implements the so-called neutral boundary condition for outflow + * boundaries. Note that this class can be used only with the stabilized + * vertex-centered discretization of the Navier-Stokes equations. + * + * This boundary condition imposes the outflow Q at boundary \f$ F \f$: + *
    + *
  • \f$ \int_Q = c \vec{u} \cdot \vec{ds} \f$, + *
+ * where \f$ \vec{ds} \f$ is the outer normal at \f$ F \f$, and + * \f$ c \f$ the concentration. + */ +template< typename TDomain> +class ConvectionDiffusionOutflowFV1 + : public ConvectionDiffusionOutflowBase +{ + private: + /// Base class type + typedef ConvectionDiffusionOutflowBase base_type; + + /// own type + typedef ConvectionDiffusionOutflowFV1 this_type; + + public: + /// World dimension + static const int dim = base_type::dim; + + public: + /// Constructor (setting default values) + ConvectionDiffusionOutflowFV1(SmartPtr< ConvectionDiffusionBase > spMaster); + + protected: + /// sets the velocity + virtual void set_velocity(SmartPtr, dim> > data) + {m_imVelocity.set_data(data);} + + + public: + /// type of trial space for each function used + virtual void prepare_setting(const std::vector& vLfeID, bool bNonRegularGrid); + + public: + /// prepares the element loop + template + void prep_elem_loop(const ReferenceObjectID roid, const int si); + + /// prepares the element for evaluation + template + void prep_elem(const LocalVector& u, GridObject* elem, const ReferenceObjectID roid, const MathVector vCornerCoords[]); + + /// finishes the element loop + template + void fsh_elem_loop(); + + /// adds the stiffness part to the local jacobian + template + void add_jac_A_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]); + + /// adds the stiffness part to the local defect + template + void add_def_A_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]); + + public: + /// dummy implementations + /// \{ + template + void add_jac_M_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]){} + template + void add_def_M_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]){} + template + void add_rhs_elem(LocalVector& d, GridObject* elem, const MathVector vCornerCoords[]){} + /// \} + + + private: + + /// adds the convective part of the local Jacobian of the momentum equation + template + inline void convective_flux_Jac + ( + const size_t ip, + const BF& bf, + LocalMatrix& J, + const LocalVector& u, + const number& stdValue + ); + /// adds the convective part of the local defect of the momentum equation + template + inline void convective_flux_defect + ( + const size_t ip, + const BF& bf, + LocalVector& d, + const LocalVector& u, + const number& stdValue + ); + + + + + protected: + /// abbreviation for pressure + static const size_t _C_ = 0; + + using base_type::m_spMaster; + using base_type::m_vBndSubSetIndex; + + /// Data import for the Velocity field + DataImport, dim > m_imVelocity; + + /// Boundary integration points of the velocity + std::vector > m_vLocIP; + std::vector > m_vGloIP; + + protected: + void register_all_funcs(bool bHang); + template + void register_func(); + +}; + +/// @} + +} // namespace ConvectionDiffusionPlugin +} // end namespace ug + +#endif /*__H__UG__PLUGINS__CONVECTION_DIFFUSION__CONVECTION_DIFFUSION_FV1__BND__OUTFLOW_FV1__*/ diff --git a/fv1/convection_diffusion_fv1.h b/fv1/convection_diffusion_fv1.h index 2891e47..68139e5 100644 --- a/fv1/convection_diffusion_fv1.h +++ b/fv1/convection_diffusion_fv1.h @@ -96,6 +96,9 @@ class ConvectionDiffusionFV1 : public ConvectionDiffusionBase /// returns the 'condensed scvf ip' flag bool condensed_FV() {return m_bCondensedFV;} + + /// returns velocity + SmartPtr, dim> > velocity() {return m_imVelocity.user_data ();} private: /// prepares assembling diff --git a/fvcr/convection_diffusion_fvcr.h b/fvcr/convection_diffusion_fvcr.h index 986d7c8..808d2ec 100644 --- a/fvcr/convection_diffusion_fvcr.h +++ b/fvcr/convection_diffusion_fvcr.h @@ -78,6 +78,9 @@ class ConvectionDiffusionFVCR : public ConvectionDiffusionBase * \param shapes upwind method */ void set_upwind(SmartPtr > shapes); + + /// returns velocity + SmartPtr, dim> > velocity() {return m_imVelocity.user_data ();} private: /// prepares the loop over all elements