diff --git a/CMakeLists.txt b/CMakeLists.txt index ca3f0b1..00fa761 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -38,6 +38,7 @@ set(pluginName ConvectionDiffusion) set(SOURCES convection_diffusion_base.cpp fv1/convection_diffusion_fv1.cpp + fv1_cutElem/convection_diffusion_fv1_cutElem.cpp fe/convection_diffusion_fe.cpp fe/convection_diffusion_stab_fe.cpp fvcr/convection_diffusion_fvcr.cpp diff --git a/convection_diffusion_plugin.cpp b/convection_diffusion_plugin.cpp index 143e155..5526961 100644 --- a/convection_diffusion_plugin.cpp +++ b/convection_diffusion_plugin.cpp @@ -36,13 +36,17 @@ #include "bridge/util.h" #include "bridge/util_domain_dependent.h" +#include "bridge/util_domain_algebra_dependent.h" #include "convection_diffusion_base.h" #include "convection_diffusion_sss.h" #include "fv1/convection_diffusion_fv1.h" +#include "fv1_cutElem/convection_diffusion_fv1_cutElem.h" #include "fe/convection_diffusion_fe.h" #include "fe/convection_diffusion_stab_fe.h" #include "fvcr/convection_diffusion_fvcr.h" #include "fv/convection_diffusion_fv.h" +#include "fv1_cutElem/diffusion_interface/diffusion_interface.h" +#include "lib_disc/spatial_disc/elem_disc/sss.h" #include "fractfv1/convection_diffusion_fractfv1.h" using namespace std; @@ -190,7 +194,20 @@ static void Domain(Registry& reg, string grp) .set_construct_as_smart_pointer(true); reg.add_class_to_group(name, "ConvectionDiffusionFV1", tag); } - + +// Convection Diffusion FV1 immersed Boundary + { + typedef ConvectionDiffusionFV1_cutElem T; + typedef ConvectionDiffusionBase TBase; + string name = string("ConvectionDiffusionFV1_cutElem").append(suffix); + reg.add_class_(name, grp) + .template add_constructor("Function(s)#Subset(s)") + .add_method("set_upwind", &T::set_upwind) + .add_method("set_testCase", &T::set_testCase) + .set_construct_as_smart_pointer(true); + reg.add_class_to_group(name, "ConvectionDiffusionFV1_cutElem", tag); + } + // Convection Diffusion FE { typedef ConvectionDiffusionFE T; @@ -243,13 +260,63 @@ static void Domain(Registry& reg, string grp) } } + /** + * Function called for the registration of Domain and Algebra dependent parts + * of the plugin. All Functions and Classes depending on both Domain and Algebra + * are to be placed here when registering. The method is called for all + * available Domain and Algebra types, based on the current build options. + * + * @param reg registry + * @param parentGroup group for sorting of functionality + */ + template + static void DomainAlgebra(Registry& reg, string grp) + { + string suffix = GetDomainAlgebraSuffix(); + string tag = GetDomainAlgebraTag(); + + typedef ApproximationSpace approximation_space_type; + typedef GridFunction function_type; + + + // ImmersedInterfaceDiffusion + { + + typedef ImmersedInterfaceDiffusion T; + typedef IImmersedInterface TBase; + string name = string("ImmersedInterfaceDiffusion").append(suffix); + reg.add_class_(name, grp) + .template add_constructor > ass, + SmartPtr > spMaster, + SmartPtr > interfaceProvider, + SmartPtr > cutElementHandler)>("domain disc, global handler") + .add_method("init", &T::init) + .add_method("set_source_data_lua", &T::set_source_data_lua) + .add_method("set_jump_data_lua", &T::set_jump_data_lua) + .add_method("set_diffusion_data_lua", &T::set_diffusion_data_lua) + .add_method("set_jump_grad_data_lua", &T::set_jump_grad_data_lua) + .add_method("get_L2Error", &T::get_L2Error) + .add_method("get_numDoFs", &T::get_numDoFs) + .add_method("set_Nitsche", &T::set_Nitsche) + .add_method("set_print_cutElemData", &T::set_print_cutElemData) + .add_method("get_numCutElements", &T::get_numCutElements) + .add_method("adjust_for_error", &T::adjust_for_error) + .add_method("initialize_threshold", &T::initialize_threshold) + .add_method("set_threshold", &T::set_threshold) + .add_method("set_analytic_solution", &T::set_analytic_solution) + .set_construct_as_smart_pointer(true); + reg.add_class_to_group(name, "ImmersedInterfaceDiffusion", tag); + } + } + + template static void Dimension(Registry& reg, string grp) { string dimSuffix = GetDimensionSuffix(); string dimTag = GetDimensionTag(); - - // singular sources and sinks + + // singular sources and sinks { typedef CDSingularSourcesAndSinks T; typedef typename T::point_sss_type TPointSSS; @@ -337,6 +404,30 @@ static void Domain(Registry& reg, string grp) .set_construct_as_smart_pointer(true); reg.add_class_to_group(name, "ConvectionDiffusionFractFV1", tag); } + + // CutElementHandler_TwoSided + { + typedef CutElementHandler_TwoSided T; + string name = string("CutElementHandler_TwoSided").append(suffix); + reg.add_class_(name, grp) + .template add_constructor mg, const char*, SmartPtr >)>("multigrid, fct names") + .set_construct_as_smart_pointer(true); + reg.add_class_to_group(name, "CutElementHandler_TwoSided", tag); + } + + // DiffusionInterfaceProvider + { + typedef DiffusionInterfaceProvider T; + string name = string("DiffusionInterfaceProvider").append(suffix); + reg.add_class_(name, grp) + .template add_constructor("") + .add_method("print", &T::print) + .add_method("add", &T::add) + .set_construct_as_smart_pointer(true); + reg.add_class_to_group(name, "DiffusionInterfaceProvider", tag); + } + + } }; // end Functionality2d3d @@ -360,6 +451,7 @@ InitUGPlugin_ConvectionDiffusion(Registry* reg, string grp) try{ RegisterDimensionDependent(*reg,grp); RegisterDomainDependent(*reg,grp); + RegisterDomainAlgebraDependent(*reg,grp); RegisterDomain2d3dDependent(*reg,grp); } UG_REGISTRY_CATCH_THROW(grp); diff --git a/fv1_cutElem/Info File 1 b/fv1_cutElem/Info File 1 new file mode 100644 index 0000000..702eb21 --- /dev/null +++ b/fv1_cutElem/Info File 1 @@ -0,0 +1,31 @@ +/* Info File for 'ImmersedInterfaceDiffusion' class: */ + + +The 'ImmersedInterfaceDiffusion' class enables the computation of elliptic problems with discontinuous coefficients, whose line along the jump forms the immersed interface. A detailed description of the mathematical theory and numerical tests can be found in [1]. + +The 'ImmersedInterfaceDiffusion' class is a derivation of the 'ImmersedInterface' class. It therefore contains all the components of an 'ImmersedInterface' class implemented according to the necessities for an immersed interface with jumping coefficients. Those are explained in more detail in the according class types. + + +The special requirements for the elliptic immersed interface implementation are the following: + +(1) The domain on BOTH sides of the interface is relevant for the physics of the problem. Therefore, the cut element computations performed by the according classes ( CutElementHandler, LocalInterfaceHandler) will be called TWICE for each original element. Accordingly, all the 'ElemDisc' methods for the computation of the local stiffnes matrix and defect contain two calls of the same methods, but get as additional parameter the inversed orientation of the interface (see convection_diffusion_fv1_cutElem.cpp). + +(2) Due to (1) two local stiffness matrices and defects will be computed during the element-local assembling process. These can NOT be stored in the usual data structure. Therefore, according local data structures are contained in the 'InterfaceHandlerLocal'. + +(3) The ansatz space used for this application is an interface adapted space. It falls into case 2 of a local-to-global mapping type (for detailed explanation, see the description in 'Info File 2' within this folder, and also 'Info File 1 ' and 'Info File 2' in the ugbase/lib_disc/spatial_disc/immersed_util folder. ). + +(4) Due to (3) additional DoFs need to be provided by the algebra for the increased, interface adapted mesh. This will be done by increasing the matrix and vector data during the initialisation pahse. + + + + +The two central methods called during the initialisation via the 'init()' method are: + +---> update_interface_data(): called by the init() of the associated CutElementHandler; computes all cut elements; +---> call of 'u.resize()'; this yields the increase of DoFs; + + + +[1] Hoellbacher S., Wittum G.: + "A sharp interface method using enriched finite elements for elliptic interface problems." + submitted to Numerische Mathematik diff --git a/fv1_cutElem/Info File 2 b/fv1_cutElem/Info File 2 new file mode 100644 index 0000000..f53b275 --- /dev/null +++ b/fv1_cutElem/Info File 2 @@ -0,0 +1,39 @@ +/* Info File for 'DiffusionInterfaceMapper' class: */ + + +Influence of the local-to-global mapping on the mathematical ansatz space: + + +On a cut element, local couplings will be computed for all the intersection points of the interface with the edges of the cut element. The mapping of these couplings to its associated global indices finally defines the property of the underlying finite element space: + +Case 1: If each coupling entry is mapped to a distinguished global index, i.e. DoF, the resulting space will be the usual finite element space on that interface adaped grid containing additional nodes along the interface and consisting (virtually) of cut elements and original elements. + +Case 2: If the coupling entry of an intersection point is mapped onto the global index of the original node, which lies accross the interface, but on the corner of the element, the finite element space will be potentially smaller, since two intersection points can share the same node accross the interface. The resulting space is a so called "flat top" space, since it results in piecewise constant solutions along the connecting line between these intersection points. + + +The 'ParticleMapper' class is of type case 2: +It maps all couplings of intersection points to the same translational (and in analogy to the rotational) global index of a particle. Therefore, the 'ParticleMapper' is a flat-top mapper. + +In contrast, the 'DiffusionInterfaceMapper' is of type case 1: +The resulting space is an interface adapted space. + + +In particular, the 'DiffusionInterfaceMapper' is a 'TWO-SIDES' mapper: As described in the 'Info File 1' within this folder, the local couplings of BOTH parts of the cut element will be mapped onto the global matrix and vectors. As a consequence, the mapping will be performed TWICE. + + +The important methods are: + +---> 'add_local_mat_to_global_interface()', which performs the local-to-global mapping first for a triangle and second for a quadrilateral, as the two parts of a 2d cut element. + +---> 'add_local_vec_to_global_interface()', which performs the local-to-global mapping first for a triangle and second for a quadrilateral, as the two parts of a 2d cut element. + + +---> modify_GlobalSol(): writes the solution within the additional interface nodes into data structures of the 'InterfaceHandlerLocal' 'in order provide global access + + + +Remark: +The 'modify_LocalData()' will NOT be used. It enables the resizing of the local data (LocalMatrix or LocalVecor) BEFORE starting the assembling on (potentially more nodes containing) cut the element. + + + diff --git a/fv1_cutElem/convection_diffusion_fv1_cutElem.cpp b/fv1_cutElem/convection_diffusion_fv1_cutElem.cpp new file mode 100644 index 0000000..5836799 --- /dev/null +++ b/fv1_cutElem/convection_diffusion_fv1_cutElem.cpp @@ -0,0 +1,1492 @@ +/* + * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt + * Author: 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__LIB_DISC__CONVECTION_DIFFUSION__CONVECTION_DIFFUSION_FV1_CUTELEM_IMPL +#define __H__UG__LIB_DISC__CONVECTION_DIFFUSION__CONVECTION_DIFFUSION_FV1_CUTELEM_IMPL + + +#include "convection_diffusion_fv1_cutElem.h" + +#include "lib_disc/spatial_disc/disc_util/fv1Cut_geom.h" +#include "lib_disc/spatial_disc/disc_util/conv_shape.h" + +namespace ug{ +namespace ConvectionDiffusionPlugin{ + + +DebugID DID_CONV_DIFF_FV1_CUTELEM("CONV_DIFF_FV1_CUTELEM"); + +enum testCase +{ + GANGL_CENTER = 0, + GANGL_OFF_CENTER, + FEDKIW_EX5, +}; + +//////////////////////////////////////////////////////////////////////////////// +// helper function +//////////////////////////////////////////////////////////////////////////////// + +/// helper function to prepare data for 'add_def_A_elem_local()' and 'add_jac_A_elem_local()' +template +template +void ConvectionDiffusionFV1_cutElem:: +get_local_data(TFVGeom& geo, const LocalVector& u, const LocalIndices& ind, LocalVector& locUOut, MathMatrix& diffusionOut, LocalVector& jumpOut, LocalVector& jump_gradOut, LocalVector& sourceOut) +{ + const bool bElementIsCut = geo.get_element_modus(); + + const int orientation = geo.get_orientation(); + + std::vector imSource; + if ( m_imSource.data_given() ) { + for ( size_t i = 0; i < 3; ++i ) + imSource.push_back(m_imSource[i]); + } + +// A. set data for non-cut elements + if ( !bElementIsCut ) + { + // set diffusion + diffusionOut *= geo.get_diffusion(geo.get_boolian_for_diffusion()); + + // set boundary condition values source, jump, jump_grad + jumpOut.resize(ind); jump_gradOut.resize(ind); + jumpOut = 0.0; jump_gradOut = 0.0; + geo.set_source(imSource, sourceOut, ind, 3, false); + + } +// B. set data for cut element + else + { + // set diffusion + diffusionOut *= geo.get_diffusion(); + + // set boundary condition values source, jump, jump_grad + int indexSize = 3; + if ( geo.get_roid() == ROID_QUADRILATERAL ) + indexSize = 4; + + // data for cut elements is stored in the InterfaceHandlerLocal and + // accessable by ElemDisc via the TFVGeom geometry class + geo.set_local_sol(locUOut, indexSize, u, orientation); + geo.set_jump_values(jumpOut, ind, indexSize); + geo.set_jump_grad_values(jump_gradOut, ind, indexSize); + geo.set_source(imSource, sourceOut, ind, indexSize, true); + } + + return; + +} + +//////////////////////////////////////////////////////////////////////////////// +// general +//////////////////////////////////////////////////////////////////////////////// + +template +ConvectionDiffusionFV1_cutElem:: +ConvectionDiffusionFV1_cutElem(const char* functions, const char* subsets) + : ConvectionDiffusionBase(functions,subsets), + m_spConvShape(new ConvectionShapesNoUpwind), + m_bNonRegularGrid(false), + m_testCase(0) +{ + register_all_funcs(m_bNonRegularGrid); +} + +template +void ConvectionDiffusionFV1_cutElem:: +prepare_setting(const std::vector& vLfeID, bool bNonRegularGrid) +{ +// check number + if(vLfeID.size() != 1) + UG_THROW("ConvectionDiffusion: Wrong number of functions given. " + "Need exactly "<<1); + + if(vLfeID[0].order() != 1 || vLfeID[0].type() != LFEID::LAGRANGE) + UG_THROW("ConvectionDiffusion FV Scheme only implemented for 1st order."); + + m_LFEID = vLfeID[0]; + +// remember + m_bNonRegularGrid = bNonRegularGrid; + +// update assemble functions + register_all_funcs(m_bNonRegularGrid); +} + +template +bool ConvectionDiffusionFV1_cutElem:: +use_hanging() const +{ + return true; +} + +//////////////////////////////////////////////////////////////////////////////// +// Assembling functions +//////////////////////////////////////////////////////////////////////////////// + +template +void ConvectionDiffusionFV1_cutElem:: +prep_assemble_loop() +{ + +} + +template +template +void ConvectionDiffusionFV1_cutElem:: +prep_elem_loop(const ReferenceObjectID roid, const int si) +{ + +// check, that upwind has been set + if(m_spConvShape.invalid()) + UG_THROW("ConvectionDiffusionFV1_cutElem::prep_elem_loop:" + " Upwind has not been set."); + +// set local positions +// if(!TFVGeom::usesHangingNodes) + if(!TFVGeom::usesHangingNodes && TFVGeom::staticLocalData) + { + static const int refDim = TElem::dim; +// TFVGeom& geo = GeomProvider::get(); + TFVGeom& geo = GeomProvider::get(m_LFEID, 1); + const MathVector* vSCVFip = geo.scvf_local_ips(); + const size_t numSCVFip = geo.num_scvf_ips(); + const MathVector* vSCVip = geo.scv_local_ips(); + const size_t numSCVip = geo.num_scv_ips(); + m_imDiffusion.template set_local_ips(vSCVFip,numSCVFip, false); + m_imVelocity.template set_local_ips(vSCVFip,numSCVFip, false); + m_imFlux.template set_local_ips(vSCVFip,numSCVFip, false); + m_imSource.template set_local_ips(vSCVip,numSCVip, false); + m_imVectorSource.template set_local_ips(vSCVFip,numSCVFip, false); + m_imReactionRate.template set_local_ips(vSCVip,numSCVip, false); + m_imReaction.template set_local_ips(vSCVip,numSCVip, false); + m_imReactionRateExpl.template set_local_ips(vSCVip,numSCVip, false); + m_imReactionExpl.template set_local_ips(vSCVip,numSCVip, false); + m_imSourceExpl.template set_local_ips(vSCVip,numSCVip, false); + m_imMassScale.template set_local_ips(vSCVip,numSCVip, false); + m_imMass.template set_local_ips(vSCVip,numSCVip, false); + + // init upwind for element type + if(!m_spConvShape->template set_geometry_type(geo)) + UG_THROW("ConvectionDiffusionFV1_cutElem::prep_elem_loop:" + " Cannot init upwind for element type."); + } +} + +template +template +void ConvectionDiffusionFV1_cutElem:: +fsh_elem_loop() +{} + +template +template +void ConvectionDiffusionFV1_cutElem:: +prep_elem(const LocalVector& u, GridObject* elem, const ReferenceObjectID roid, const MathVector vCornerCoords[]) +{ +// Update Geometry for this element +//static TFVGeom& geo = GeomProvider::get(); + TFVGeom& geo = GeomProvider::get(m_LFEID,1); + +// fix: set orientation initially globally! + geo.set_orientation(1); + + + try{ + UG_DLOG(DID_CONV_DIFF_FV1_CUTELEM, 2, ">>OCT_DISC_DEBUG: " << "convection_diffusion_fv1.cpp: " << "prep_elem(): update(): "<< roid << std::endl); + geo.update(elem, vCornerCoords, &(this->subset_handler())); + }UG_CATCH_THROW("ConvectionDiffusionFV1_cutElem::prep_elem:" + " Cannot update Finite Volume Geometry."); + +// set local positions +// if(TFVGeom::usesHangingNodes) + if(TFVGeom::usesHangingNodes || !TFVGeom::staticLocalData) + { + const int refDim = TElem::dim; + const MathVector* vSCVFip = geo.scvf_local_ips(); + const size_t numSCVFip = geo.num_scvf_ips(); + const MathVector* vSCVip = geo.scv_local_ips(); + const size_t numSCVip = geo.num_scv_ips(); + m_imDiffusion.template set_local_ips(vSCVFip,numSCVFip); + m_imVelocity.template set_local_ips(vSCVFip,numSCVFip); + m_imFlux.template set_local_ips(vSCVFip,numSCVFip); + m_imSource.template set_local_ips(vSCVip,numSCVip); + m_imVectorSource.template set_local_ips(vSCVFip,numSCVFip); + m_imReactionRate.template set_local_ips(vSCVip,numSCVip); + m_imReaction.template set_local_ips(vSCVip,numSCVip); + m_imReactionRateExpl.template set_local_ips(vSCVip,numSCVip); + m_imReactionExpl.template set_local_ips(vSCVip,numSCVip); + m_imSourceExpl.template set_local_ips(vSCVip,numSCVip); + m_imMassScale.template set_local_ips(vSCVip,numSCVip); + m_imMass.template set_local_ips(vSCVip,numSCVip); +/* + if(m_spConvShape.valid()) + if(!m_spConvShape->template set_geometry_type(geo)) + UG_THROW("ConvectionDiffusionFV1_cutElem::prep_elem_loop:" + " Cannot init upwind for element type."); + */ + } + + // set global positions + const MathVector* vSCVFip = geo.scvf_global_ips(); + const size_t numSCVFip = geo.num_scvf_ips(); + const MathVector* vSCVip = geo.scv_global_ips(); + const size_t numSCVip = geo.num_scv_ips(); + m_imDiffusion. set_global_ips(vSCVFip, numSCVFip); + m_imVelocity. set_global_ips(vSCVFip, numSCVFip); + m_imFlux. set_global_ips(vSCVFip, numSCVFip); + m_imSource. set_global_ips(vSCVip, numSCVip); + m_imVectorSource. set_global_ips(vSCVFip, numSCVFip); + m_imReactionRate. set_global_ips(vSCVip, numSCVip); + m_imReactionRateExpl. set_global_ips(vSCVip, numSCVip); + m_imReactionExpl. set_global_ips(vSCVip, numSCVip); + m_imSourceExpl. set_global_ips(vSCVip, numSCVip); + m_imReaction. set_global_ips(vSCVip, numSCVip); + m_imMassScale. set_global_ips(vSCVip, numSCVip); + m_imMass. set_global_ips(vSCVip, numSCVip); +} + +template +static TVector CalculateCenter(GridObject* o, const TVector* coords) +{ + TVector v; + VecSet(v, 0); + + size_t numCoords = 0; + switch(o->base_object_id()){ + case VERTEX: numCoords = 1; break; + case EDGE: numCoords = static_cast(o)->num_vertices(); break; + case FACE: numCoords = static_cast(o)->num_vertices(); break; + case VOLUME: numCoords = static_cast(o)->num_vertices(); break; + default: UG_THROW("Unknown element type."); break; + } + + for(size_t i = 0; i < numCoords; ++i) + VecAdd(v, v, coords[i]); + + if(numCoords > 0) + VecScale(v, v, 1. / (number)numCoords); + + return v; +} + +template +template +void ConvectionDiffusionFV1_cutElem:: +add_jac_A_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]) +{ +// get finite volume geometry + TFVGeom& geo = GeomProvider::get(m_LFEID,1); + +// get the modus of the current element + const bool bElementIsCut = geo.get_element_modus(); + +//////////////////////////////////////////////////////////////////////////// +// (A) usual assembling for non-cut element: + if ( !bElementIsCut ) + { + LocalVector dummyU = u; + this->template add_jac_A_elem_local (geo, u, J, dummyU, elem, vCornerCoords); + return; + } + +//////////////////////////////////////////////////////////////////////////// +// (B) call 'add_jac_A_elem_local_cut()', which handles the adaptions +// necessary for the cut element + +// initialize cut element data: + geo.L2Error_init(); + geo.resize_local_data(u); + + +//////////////////////////////////////////////////////////////////////////// +// (B1) First: call of 'add_jac_A_elem_local_cut()' with orientation = 1: + int orientation = 1; + geo.set_orientation(orientation); + + add_jac_A_elem_local_cut(geo, u, elem, vCornerCoords); + +//////////////////////////////////////////////////////////////////////////// +// (B2) Second: call of 'add_jac_A_elem_local_cut()' with orientation = -1: + orientation *= -1; + geo.set_orientation(orientation); + +// recompute local data on cut element due to new orientation: the local +// TFVGeom geometry data of the other part of the cut element will be computed + try{ + geo.update(elem, vCornerCoords, &(this->subset_handler())); + }UG_CATCH_THROW("ConvectionDiffusionFV1_cutElem::update:" + " Cannot update Finite Volume Geometry."); + + add_jac_A_elem_local_cut(geo, u, elem, vCornerCoords); + + +} + +template +template +void ConvectionDiffusionFV1_cutElem:: +add_jac_A_elem_local_cut(TFVGeom& geo, const LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]) +{ + bool boundary = false; + +// shiftTag = true: in case of double DoFs on interface (for strong discontinuity on the interface) + bool shiftTag = geo.get_bScaleDoFs(); + + ReferenceObjectID roidCheck = geo.get_roid(); + int indexSize = 3; + if ( roidCheck == ROID_QUADRILATERAL ) indexSize = 4; + +// get the pointer for access and storage into 'InterfaceHandlerLocalDiffusion' class + LocalMatrix& locJ = geo.get_jacobian(roidCheck); + LocalVector& locU = geo.get_solution(roidCheck); + + locJ = 0; + +// finally: compute the local stiffness matrix on the cut element + this->template add_jac_A_elem_local (geo, u, locJ, locU, elem, vCornerCoords); + +// compute the boundary condition on the interface of the cut element + if ( boundary ) + { + geo.reset_jacobian_on_interface(locJ, indexSize); + this->template add_jac_A_elem_boundary (geo, locJ, locU, elem, vCornerCoords); + } + +// write computed local stiffness matrix to data storage in class 'InterfaceHandlerLocalDiffusion' + geo.set_jacobian(locJ, roidCheck); + +// shiftTag necessary for local to global mapping + geo.set_DoF_tag(shiftTag, roidCheck); + + +} + + +template +template +void ConvectionDiffusionFV1_cutElem:: +add_jac_A_elem_local(TFVGeom& geo, const LocalVector& u, LocalMatrix& J, LocalVector& locU, + GridObject* elem, const MathVector vCornerCoords[]) +{ +//////////////////////////////////////////////////////////////////////// +// Pre-processing for assembling on cut Element: + + MathMatrix diffusion; + MatSet(diffusion, 0); + MatDiagSet(diffusion, 1.0); + + LocalVector source, jump; + +// get the local data (jump, source) on the cut element + get_local_data(geo, u, u.get_indices(), locU, diffusion, jump, jump, jump); + +//////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////// +// NOW: Standard ConvectionDiffusionFV1-assembling starts here: + + +// Diff. Tensor times Gradient + MathVector Dgrad; + +// get conv shapes + const IConvectionShapes& convShape = get_updated_conv_shapes(geo); + +// Diffusion and Velocity Term + if(m_imDiffusion.data_given() || m_imVelocity.data_given()) + { + // loop Sub Control Volume Faces (SCVF) + for(size_t ip = 0; ip < geo.num_scvf(); ++ip) + { + // get current SCVF + const typename TFVGeom::SCVF& scvf = geo.scvf(ip); + + //////////////////////////////////////////////////// + // Diffusive Term + //////////////////////////////////////////////////// + if(m_imDiffusion.data_given()) + { + // DID_CONV_DIFF_FV1_CUTELEM + number D_diff_flux_sum = 0.0; + + // loop shape functions + for(size_t sh = 0; sh < scvf.num_sh(); ++sh) + { + // Compute Diffusion Tensor times Gradient + MatVecMult(Dgrad, diffusion, scvf.global_grad(sh)); + + // Compute flux at IP + const number D_diff_flux = VecDot(Dgrad, scvf.normal()); + UG_DLOG(DID_CONV_DIFF_FV1_CUTELEM, 2, ">>OCT_DISC_DEBUG: " << "convection_diffusion_fv1.cpp: " << "add_jac_A_elem(): " << "sh # " << sh << " ; normalSize scvf # " << ip << ": " << VecLength(scvf.normal()) << "; \t from "<< scvf.from() << "; to " << scvf.to() << "; D_diff_flux: " << D_diff_flux << "; scvf.global_grad(sh): " << scvf.global_grad(sh) << std::endl); + + // Add flux term to local matrix // HIER MATRIXINDIZES!!! + UG_ASSERT((scvf.from() < J.num_row_dof(_C_)) && (scvf.to() < J.num_col_dof(_C_)), + "Bad local dof-index on element with object-id " << elem->base_object_id() + << " with center: " << CalculateCenter(elem, vCornerCoords)); + + J(_C_, scvf.from(), _C_, sh) -= D_diff_flux; + J(_C_, scvf.to() , _C_, sh) += D_diff_flux; + + // DID_CONV_DIFF_FV1_CUTELEM + D_diff_flux_sum += D_diff_flux; + } + + UG_DLOG(DID_CONV_DIFF_FV1_CUTELEM, 2, "D_diff_flux_sum = " << D_diff_flux_sum << std::endl << std::endl); + } + + + ///////////////////////////////////////////////////////////////////////////// + // Additional diffusive Term due to jump in solution at the interface + // u^+ - u^- = jump + ///////////////////////////////////////////////////////////////////////////// + if ( 0 ) + { + // loop shape functions + for(size_t sh = 0; sh < scvf.num_sh(); ++sh) + { + // Compute Diffusion Tensor times Gradient + MatVecMult(Dgrad, diffusion, scvf.global_grad(sh)); + + // Compute flux at IP + const number D_diff_flux = jump(_C_,sh) * VecDot(Dgrad, scvf.normal()); + + J(_C_, scvf.from(), _C_, sh) -= D_diff_flux; + J(_C_, scvf.to() , _C_, sh) += D_diff_flux; + } + + } + + //////////////////////////////////////////////////// + // Convective Term + //////////////////////////////////////////////////// + if(m_imVelocity.data_given()) + { + // Add Flux contribution + for(size_t sh = 0; sh < convShape.num_sh(); ++sh) + { + const number D_conv_flux = convShape(ip, sh); + + // Add flux term to local matrix + J(_C_, scvf.from(), _C_, sh) += D_conv_flux; + J(_C_, scvf.to(), _C_, sh) -= D_conv_flux; + } + } + + // no explicit dependency on flux import + } + } + + //UG_LOG("Local Matrix is: \n"< +template +void ConvectionDiffusionFV1_cutElem:: +add_jac_A_elem_boundary(TFVGeom& geo, LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]) +{ +// get parameter + double diffusion = 0.0; + if ( !geo.get_element_modus() ) + diffusion = geo.get_diffusion(geo.get_boolian_for_diffusion()); + else + diffusion = geo.get_diffusion(); + +// get inner boundary faces for assembling (stored in 'InterfaceHandlerLocalDiffusion' class) + std::vector& vBF = geo.get_boundary_faces(); + if ( vBF.size() > 2 ) + UG_THROW("add_jac_A_elem_boundary(): vBF.size() is greater than 2: " << vBF.size() << "\n"); + + if ( vBF.size() == 0 && geo.get_element_modus() ) + UG_THROW("add_jac_A_elem_boundary(): Error! If vBF.size() == 0, then the element is cut!\n"); + +//////////////////////////////////////////////////////////////////////////////// +// REMARK: NO loop integration points! +// /* for(size_t ip = 0; ip < vBF.size(); ++ip) */ +// Reason: the length of the normal is already the length of the total face (NOT the scvf!) +//////////////////////////////////////////////////////////////////////////////// + +// loop integration points + for(size_t ip = 0; ip < vBF.size(); ++ip) + { + typename TFVGeom::BF bf = vBF[ip]; + + // loop trial space + for(size_t sh = 0; sh < bf.num_sh(); ++sh) + { + // add to local matrix + J(_C_, bf.node_id(), _C_, sh) += VecDot(bf.global_grad(sh), bf.normal()); + J(_C_, bf.node_id(), _C_, sh) *= diffusion; + } + } + +} + + +template +template +void ConvectionDiffusionFV1_cutElem:: +add_jac_M_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]) +{ +// get finite volume geometry + static const TFVGeom& geo = GeomProvider::get(m_LFEID,1); + + if(!m_imMassScale.data_given()) return; + +// loop Sub Control Volumes (SCV) + for(size_t ip = 0; ip < geo.num_scv(); ++ip) + { + // get current SCV + const typename TFVGeom::SCV& scv = geo.scv(ip); + + // get associated node + const int co = scv.node_id(); + + // Add to local matrix + J(_C_, co, _C_, co) += scv.volume() * m_imMassScale[ip]; + } + +} + +template +template +void ConvectionDiffusionFV1_cutElem:: +add_def_A_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]) +{ +// get finite volume geometry + TFVGeom& geo = GeomProvider::get(m_LFEID,1); + +// get the modus of the current element + const bool bElementIsCut = geo.get_element_modus(); + + std::vector imSource; + if ( m_imSource.data_given() ) { + for ( size_t i = 0; i < 3; ++i ) + imSource.push_back(m_imSource[i]); + } + +//////////////////////////////////////////////////////////////////////// +// (A) usual assembling for non-cut elements + if ( !bElementIsCut ) + { + LocalVector dummyU = u; + this->template add_def_A_elem_local (geo, u, d, dummyU, elem, vCornerCoords); + + number intValElem = this->template add_l2error_A_elem (geo, ROID_TRIANGLE, d, u, elem); + geo.L2Error_add(intValElem); + + return; + } + +///////////////////////////////////////////////////////////////////////////// +// (B) call 'add_def_A_elem_local_cut()', which handles the adaptions +// due to the cut element + +// initialize cut element data: + geo.resize_local_data(u); + +////////////////////////////////////////////////////////////////////////////// +// (B1) First: call of 'add_def_A_elem_local_cut()' with orientation = 1: + int orientation = 1; + geo.set_orientation(orientation); + + add_def_A_elem_local_cut (geo, u, elem, vCornerCoords); + + +///////////////////////////////////////////////////////////////////////////// +// (B2) Second: call of 'add_def_A_elem_local_cut()' with orientation = -1: + orientation *= -1; + geo.set_orientation(orientation); + +// recompute local data on cut element due to new orientation: the local +// TFVGeom geometry data of the other part of the cut element will be computed + try{ + geo.update(elem, vCornerCoords, &(this->subset_handler())); + }UG_CATCH_THROW("ConvectionDiffusionFV1_cutElem::update:" + " Cannot update Finite Volume Geometry."); + + add_def_A_elem_local_cut (geo, u, elem, vCornerCoords); + + +} + +template +template +void ConvectionDiffusionFV1_cutElem:: +add_def_A_elem_local_cut(TFVGeom& geo, const LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]) +{ + bool boundary = false; + bool add = true; + +// get data + ReferenceObjectID roidCheck = geo.get_roid(); + bool shiftTag = geo.get_bScaleDoFs(); // bScaleDoFs = true: in case of double DoFs on interface! + int indexSize = 3; + if ( roidCheck == ROID_QUADRILATERAL ) indexSize = 4; + +// get the pointer for access and storage into 'InterfaceHandlerLocalDiffusion' class + LocalVector& locD = geo.get_defect(roidCheck); + LocalVector& locU = geo.get_solution(roidCheck); + + locD = 0; + +// finally: compute the local defect on the cut element + this->template add_def_A_elem_local (geo, u, locD, locU, elem, vCornerCoords); + +// compute the boundary condition on the interface of the cut element + if ( boundary ) + { + geo.reset_defect_on_interface(locD, indexSize); + this->template add_def_A_elem_boundary (geo, locD, locU, elem, vCornerCoords); + } + +// write computed local defect to data storage in class 'InterfaceHandlerLocalDiffusion' + geo.set_defect(locD, roidCheck); + +// shiftTag necessary for local to global mapping + geo.set_DoF_tag(shiftTag, roidCheck); + +// compute the l2Error on the cut element + number intValElem = this->template add_l2error_A_elem (geo, roidCheck, locD, locU, elem); + if ( add ) geo.L2Error_add(intValElem); + +} + +template +template +void ConvectionDiffusionFV1_cutElem:: +add_def_A_elem_local(TFVGeom& geo, const LocalVector& u, LocalVector& d, LocalVector& locU, GridObject* elem, const MathVector vCornerCoords[]) +{ +//////////////////////////////////////////////////////////////////////// +// Pre-processing for assembling on cut Element: + + MathMatrix diffusion; + MatSet(diffusion, 0); + MatDiagSet(diffusion, 1.0); + + LocalVector jump, jump_grad, source; + +// get the local data (jump, jump_grad, source) on the cut element + get_local_data(geo, u, u.get_indices(), locU, diffusion, jump, jump_grad, source); + +//////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////// +// NOW: Standard ConvectionDiffusionFV1-assembling starts here: + + if(m_imDiffusion.data_given() || m_imVelocity.data_given() || m_imFlux.data_given()) + { + // loop Sub Control Volume Faces (SCVF) + for(size_t ip = 0; ip < geo.num_scvf(); ++ip) + { + // get current SCVF + const typename TFVGeom::SCVF& scvf = geo.scvf(ip); + + ///////////////////////////////////////////////////// + // Diffusive Term + ///////////////////////////////////////////////////// + if(m_imDiffusion.data_given()) + { + // to compute D \nabla c + MathVector Dgrad_c, grad_c; + + // compute gradient and shape at ip + VecSet(grad_c, 0.0); + for(size_t sh = 0; sh < scvf.num_sh(); ++sh) + VecScaleAppend(grad_c, locU(_C_,sh), scvf.global_grad(sh)); + + // scale by diffusion tensor + MatVecMult(Dgrad_c, diffusion, grad_c); + + // Compute flux + const number diff_flux = VecDot(Dgrad_c, scvf.normal()); + + // Add to local defect + d(_C_, scvf.from()) -= diff_flux; + d(_C_, scvf.to() ) += diff_flux; + } + + ///////////////////////////////////////////////////////////////////////////// + // Additional diffusive Term due to jump in solution at the interface + // u^+ - u^- = jump + ///////////////////////////////////////////////////////////////////////////// + + // to compute D \nabla c=Id_interface + MathVector Dgrad, grad; + + // compute gradient and shape at ip + VecSet(grad, 0.0); + for(size_t sh = 0; sh < scvf.num_sh(); ++sh) + VecScaleAppend(grad, jump(_C_,sh), scvf.global_grad(sh)); + + // scale by diffusion tensor + MatVecMult(Dgrad, diffusion, grad); + + // Compute flux + const number diff_flux = VecDot(Dgrad, scvf.normal()); + + // Add to local defect + d(_C_, scvf.from()) -= diff_flux; + d(_C_, scvf.to() ) += diff_flux; + + } + + } + + ///////////////////////////////////////////////////// + // add rhs during same method! + // --> in elem_disc_assemble_util, the method 'add_rhs_elem()' + // adds the local vector otherwise! NOT functional!! + ///////////////////////////////////////////////////// + +// loop integration points + for ( size_t ip = 0; ip < geo.num_scv(); ++ip ) + { + // get current SCV + const typename TFVGeom::SCV& scv = geo.scv( ip ); + + // get associated node + int co = scv.node_id(); + + // Add to local rhs + d(_C_, co) -= source(_C_, co) * scv.volume(); + } + + + ///////////////////////////////////////////////////////////////////////////// + // Additional source Term due to jump in gradient at the interface + // (\nabla u^+ - \nabla u^-)\cdot n = h(x) * |n| + ///////////////////////////////////////////////////////////////////////////// + + std::vector& vBF = geo.get_boundary_faces(); + +// some checks: + if ( vBF.size() > 2 ) + UG_THROW("add_def_A_elem_local(): vBF.size() is greater than 2: " << vBF.size() << "\n"); + if ( vBF.size() == 0 && geo.get_element_modus() ) + UG_THROW("add_def_A_elem_local(): Error! If vBF.size() == 0, then the element is cut!\n"); + +// loop integration points + for(size_t ip = 0; ip < vBF.size(); ++ip) + { + typename TFVGeom::BF bf = vBF[ip]; + + // Add to local rhs + d(_C_, bf.node_id()) += jump_grad(_C_,bf.node_id()) * bf.Vol; + } + + +} + +template +template +void ConvectionDiffusionFV1_cutElem:: +add_def_A_elem_boundary(TFVGeom& geo, LocalVector& d, LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]) +{ +// get parameter + double diffusion = 0.0; + if ( !geo.get_element_modus() ) + diffusion = geo.get_diffusion(geo.get_boolian_for_diffusion()); + else + diffusion = geo.get_diffusion(); + +// get inner boundary faces for assembling (stored in 'InterfaceHandlerLocalDiffusion' class) + std::vector& vBF = geo.get_boundary_faces(); + if ( vBF.size() > 2 ) + UG_THROW("add_def_A_elem(): vBF.size() is greater than 2: " << vBF.size() << "\n"); + + if ( vBF.size() == 0 && geo.get_element_modus() ) + UG_THROW("add_def_A_elem_boundary(): Error! If vBF.size() == 0, then the element is cut!\n"); + + MathVector Dgrad; + VecSet(Dgrad, 0.0); + + +//////////////////////////////////////////////////////////////////////////////// +// REMARK: NO loop integration points! +// /* for(size_t ip = 0; ip < vBF.size(); ++ip) */ +// Reason: the length of the normal is already the length of the total face (NOT the scvf!) +//////////////////////////////////////////////////////////////////////////////// + +// loop integration points + for(size_t ip = 0; ip < vBF.size(); ++ip) + { + typename TFVGeom::BF bf = vBF[ip]; + VecSet(Dgrad, 0.0); + + // loop trial space + for(size_t sh = 0; sh < bf.num_sh(); ++sh) + { + // Diffusion + VecScaleAppend(Dgrad, diffusion * u(_C_, sh), bf.global_grad(sh)); + } + + // add to local vector + d(_C_, bf.node_id()) += VecDot(Dgrad, bf.normal()); + + } + +} + + +//////////////////////////////////////////////////////////////////////////////// +/// +/// methods for cut element error computation via 'add_l2error_A_elem()' +/// --> hard coded c++ functions for some numerical test examples +/// (not optimal handling of user data) +/// --> the choice of different test examples can be handled via lua: +/// --> lua-call: 'elemDisc:set_testCase(int testCase)' +/// testCase-enumerator, see top of that file! +/// +//////////////////////////////////////////////////////////////////////////////// + +template +number get_exact_sol_test(MathVector position) +{ + return sin(2*M_PI*position[0]) + sin(2*M_PI*position[1]); +} + + + +template +number get_exact_sol_Gangl(MathVector position, size_t testCase) +{ + MathVector<2> center(0.0); + if (testCase == GANGL_OFF_CENTER ) + { + center[0] = -0.08; + center[1] = 0.3; + } + + double kappa_2 = 10.0; + double dist_x = position[0] - center[0]; + double dist_y = position[1] - center[1]; + double sqR = 0.4*0.4; + + double sqDist = dist_x*dist_x + dist_y*dist_y; + double dist = sqrt(sqDist); + + double returnValue = -4*kappa_2*kappa_2*sqR*sqDist + 2*sqR*sqR*kappa_2*(2*kappa_2 - 1); + + if ( dist >= 0.4 ) + returnValue = -2*kappa_2*sqDist*sqDist; + + return returnValue; +} + +template +MathVector get_exact_grad_Gangl(MathVector position, size_t testCase) +{ + MathVector<2> center(0.0); + if (testCase == GANGL_OFF_CENTER ) + { + center[0] = -0.08; + center[1] = 0.3; + } + + double kappa_2 = 10.0; + double dist_x = position[0] - center[0]; + double dist_y = position[1] - center[1]; + double sqR = 0.4*0.4; + + double sqDist = dist_x*dist_x + dist_y*dist_y; + double dist = sqrt(sqDist); + + double factor = -8*kappa_2*kappa_2*sqR; + + if ( dist >= 0.4 ) + factor = -8*kappa_2*sqDist; + + MathVector returnVector; + returnVector[0] = factor*dist_x; + returnVector[1] = factor*dist_y; + + return returnVector; +} + +template +MathVector get_exact_grad_FedkiwEx5(MathVector position) +{ + double center_x = 0.0; + double center_y = 0.0; + double radius = 0.5; + + double dist_x = position[0] - center_x; + double dist_y = position[1] - center_y; + + double sqDist = dist_x*dist_x + dist_y*dist_y; + double dist = sqrt(sqDist); + + double absValue = position[0]*position[0] + position[1]* position[1]; + + MathVector returnVector; + returnVector[0] = 0.0; + returnVector[1] = 0.0; + + double factor = 1.0/absValue; + if ( dist >= radius ) + { + returnVector[0] = factor*position[0]; + returnVector[1] = factor*position[1]; + } + + return returnVector; + +} + +template +number get_exact_sol_FedkiwEx5(MathVector position) +{ + double center_x = 0.0; + double center_y = 0.0; + + double dist_x = position[0]-center_x; + double dist_y = position[1]-center_y; + + double sqDist = dist_x*dist_x + dist_y*dist_y; + double dist = sqrt(sqDist); + + double returnValue = 1.0; + + if ( dist > 0.5 ) + returnValue = 1.0 + log(2*dist); + + return returnValue; +} + +template +number get_exact_sol_FedkiwEx3(MathVector position) +{ + double center_x = 0.5; + double center_y = 0.5; + double radius = 0.25; + + double dist_x = position[0] - center_x; + double dist_y = position[1] - center_y; + + double sqDist = dist_x*dist_x + dist_y*dist_y; + double dist = sqrt(sqDist); + + double absValue = position[0]*position[0] + position[1]* position[1]; + + double returnValue = 0.0; + + if ( dist <= radius ) + returnValue = exp(-absValue); + + return returnValue; +} + + +////////////////////////////////////////////////////////////////////// +// code see ugbase/lib_disc/function_spaces/integrate.h: evaluate() for +// --> L2ErrorIntegrand (for value) +// --> H1ErrorIntegrand (for gradient): lines 1873-1910 +// +// called bei Integrate() via method 'integrand.values': +// integrand.values(&(vValue[0]), &(vGlobIP[0]), +// pElem, &vCorner[0], rQuadRule.points(), +// &(vJT[0]), +// numIP); +// +////////////////////////////////////////////////////////////////////// +template +template +number ConvectionDiffusionFV1_cutElem:: +add_l2error_A_elem(TFVGeom& geo, ReferenceObjectID roid, LocalVector& d, const LocalVector& u, GridObject* elem) +{ + // get the flag for the choice of numerical example to be used as analytical solution + size_t testCase = get_testCase(); + + std::vector > vCorner; + std::vector > vGlobIP; + std::vector > vLocIP; + std::vector > vJT; + std::vector vValue; + std::vector vValueGrad; + + QuadType type = GetQuadratureType("best"); + + const QuadratureRule& rQuadRule + = QuadratureRuleProvider::get(roid, 1, type); + + // get reference element mapping by reference object id + DimReferenceMapping& mapping + = ReferenceMappingProvider::get(roid); + + // number of integration points + const size_t numIP = rQuadRule.size(); + + // get all corner coordinates + // CollectCornerCoordinates(vCorner, *pElem, aaPos, true); + + const DimReferenceElement& rRefElem + = ReferenceElementProvider::get(roid); + + vCorner.clear(); + // remember global position of nodes + for(size_t i = 0; i < rRefElem.num(0); ++i) + vCorner.push_back(geo.get_corner(i)); + + + // update the reference mapping for the corners + mapping.update(vCorner); + + // compute global integration points + vGlobIP.resize(numIP); + mapping.local_to_global(&(vGlobIP[0]), rQuadRule.points(), numIP); + + // compute local integration points + vLocIP.resize(numIP); + for(size_t ip = 0; ip < numIP; ++ip) + vLocIP[ip] = rQuadRule.points()[ip]; + + + // compute transformation matrices + vJT.resize(numIP); + mapping.jacobian_transposed(&(vJT[0]), rQuadRule.points(), numIP); + + const size_t num_sh = geo.num_scvf(); + + if ( num_sh != rRefElem.num(0) ) + UG_THROW("wrong number of corners: sh = " << num_sh << ", but rRefElem.num(0) = " << rRefElem.num(0) << "\n"); + + // compute integrand values at integration points + vValue.resize(numIP); + vValueGrad.resize(numIP); + + number exactSolIP; + MathVector exactGradIP; + try + { + // loop all integration points + for(size_t ip = 0; ip < numIP; ++ip) + { + // compute exact solution at integration point + if ( testCase == GANGL_CENTER || testCase == GANGL_OFF_CENTER ) + exactSolIP = get_exact_sol_Gangl(vGlobIP[ip], testCase); + else if ( testCase == FEDKIW_EX5 ) + exactSolIP = get_exact_sol_FedkiwEx5(vGlobIP[ip]); + + + // compute exact gradient at integration point + if ( testCase == GANGL_CENTER || testCase == GANGL_OFF_CENTER ) + exactGradIP = get_exact_grad_Gangl(vGlobIP[ip], testCase); + else if ( testCase == FEDKIW_EX5 ) + exactGradIP = get_exact_grad_FedkiwEx5(vGlobIP[ip]); + + + // compute approximated solution at integration point + number approxSolIP = 0.0; + MathVector locTmp; VecSet(locTmp, 0.0); + + const typename TFVGeom::SCV& scv = geo.scv(ip); + + for(size_t sh = 0; sh < num_sh; ++sh) + { + // add shape fct at ip * value at shape + approxSolIP += u(_C_,sh) * geo.get_shape(sh, vLocIP[ip], roid); + + // add gradient at ip + VecScaleAppend(locTmp, u(_C_,sh), scv.local_grad(sh)); + } + + // get squared of difference + vValue[ip] = (exactSolIP - approxSolIP); + vValue[ip] *= vValue[ip]; + + // compute global gradient + MathVector approxGradIP; + MathMatrix JTInv; + Inverse(JTInv, vJT[ip]); + MatVecMult(approxGradIP, JTInv, locTmp); + + // get error of gradient + vValueGrad[ip] = VecDistanceSq(approxGradIP, exactGradIP); + + + } + /* integrand.values(&(vValue[0]), &(vGlobIP[0]), + pElem, &vCorner[0], rQuadRule.points(), + &(vJT[0]), + numIP); + */ + } + UG_CATCH_THROW("Unable to compute values of integrand at integration point."); + + // reset contribution of this element + number intValElem = 0; + + // loop integration points + for(size_t ip = 0; ip < numIP; ++ip) + { + // get quadrature weight + const number weightIP = rQuadRule.weight(ip); + + // get determinate of mapping + const number det = SqrtGramDeterminant(vJT[ip]); + + // add contribution of integration point + intValElem += vValue[ip] * weightIP * det; + // intValElem += vValueGrad[ip] * weightIP * det; + + } + + // add to global sum + + return intValElem; +} + + +template +template +void ConvectionDiffusionFV1_cutElem:: +add_def_A_expl_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]) +{ +// get finite volume geometry +// static const TFVGeom& geo = GeomProvider::get(); + static const TFVGeom& geo = GeomProvider::get(m_LFEID,1); + +// reaction rate + if(m_imReactionRateExpl.data_given()) + { + // loop Sub Control Volumes (SCV) + for(size_t ip = 0; ip < geo.num_scv(); ++ip) + { + // get current SCV + const typename TFVGeom::SCV& scv = geo.scv(ip); + + // get associated node + const int co = scv.node_id(); + + // Add to local defect + d(_C_, co) += u(_C_, co) * m_imReactionRateExpl[ip] * scv.volume(); + } + } + +// reaction + if(m_imReactionExpl.data_given()) + { + // loop Sub Control Volumes (SCV) + for(size_t ip = 0; ip < geo.num_scv(); ++ip) + { + // get current SCV + const typename TFVGeom::SCV& scv = geo.scv(ip); + + // get associated node + const int co = scv.node_id(); + + // Add to local defect + d(_C_, co) += m_imReactionExpl[ip] * scv.volume(); + } + } + + if(m_imSourceExpl.data_given()) + { + // loop Sub Control Volumes (SCV) + for(size_t ip = 0; ip < geo.num_scv(); ++ip) + { + // get current SCV + const typename TFVGeom::SCV& scv = geo.scv(ip); + + // get associated node + const int co = scv.node_id(); + + // Add to local rhs + d(_C_, co) -= m_imSourceExpl[ip] * scv.volume(); + } + } +} + +template +template +void ConvectionDiffusionFV1_cutElem:: +add_def_M_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]) +{ +// get finite volume geometry +// static const TFVGeom& geo = GeomProvider::get(); + static const TFVGeom& geo = GeomProvider::get(m_LFEID,1); + + if(!m_imMassScale.data_given() && !m_imMass.data_given()) return; + +// loop Sub Control Volumes (SCV) + for(size_t ip = 0; ip < geo.num_scv(); ++ip) + { + // get current SCV + const typename TFVGeom::SCV& scv = geo.scv(ip); + + // get associated node + const int co = scv.node_id(); + + // mass value + number val = 0.0; + + // multiply by scaling + if(m_imMassScale.data_given()) + val += m_imMassScale[ip] * u(_C_, co); + + // add mass + if(m_imMass.data_given()) + val += m_imMass[ip]; + + // Add to local defect + d(_C_, co) += val * scv.volume(); + } +} + + +template +template +void ConvectionDiffusionFV1_cutElem:: +add_rhs_elem(LocalVector& d, GridObject* elem, const MathVector vCornerCoords[]) +{ + ///////////////////////////////////////////////////// + // add rhs ALLREADY(!) during 'add_def_A_elem_local()' method! + // --> in elem_disc_assemble_util, the method 'add_rhs_elem()' adds the local vector otherwise! NOT functional!! + ///////////////////////////////////////////////////// + + return; + + // get finite volume geometry + // static const TFVGeom& geo = GeomProvider::get(); + static const TFVGeom& geo = GeomProvider::get(m_LFEID,1); + + // loop Sub Control Volumes (SCV) + if ( m_imSource.data_given() ) { + for ( size_t ip = 0; ip < geo.num_scv(); ++ip ) { + // get current SCV + const typename TFVGeom::SCV& scv = geo.scv( ip ); + + // get associated node + const int co = scv.node_id(); + + // Add to local rhs + d(_C_, co) += m_imSource[ip] * scv.volume(); + } + } + + // loop Sub Control Volumes (SCVF) + if ( m_imVectorSource.data_given() ) { + for ( size_t ip = 0; ip < geo.num_scvf(); ++ip ) { + // get current SCVF + const typename TFVGeom::SCVF& scvf = geo.scvf( ip ); + + // Add to local rhs + d(_C_, scvf.from()) -= VecDot(m_imVectorSource[ip], scvf.normal() ); + d(_C_, scvf.to() ) += VecDot(m_imVectorSource[ip], scvf.normal() ); + } + } +} + + + +//////////////////////////////////////////////////////////////////////////////// +// upwind +//////////////////////////////////////////////////////////////////////////////// + +template +void ConvectionDiffusionFV1_cutElem:: +set_upwind(SmartPtr > shapes) {m_spConvShape = shapes;} + +// computes the linearized defect w.r.t to the velocity +template +const typename ConvectionDiffusionFV1_cutElem::conv_shape_type& +ConvectionDiffusionFV1_cutElem:: +get_updated_conv_shapes(const FVGeometryBase& geo) +{ +// compute upwind shapes for transport equation +// \todo: we should move this computation into the preparation part of the +// disc, to only compute the shapes once, reusing them several times. + if(m_imVelocity.data_given()) + { + // get diffusion at ips + const MathMatrix* vDiffusion = NULL; + if(m_imDiffusion.data_given()) vDiffusion = m_imDiffusion.values(); + + // update convection shapes + if(!m_spConvShape->update(&geo, m_imVelocity.values(), vDiffusion, true)) + { + UG_LOG("ERROR in 'ConvectionDiffusionFV1_cutElem::get_updated_conv_shapes': " + "Cannot compute convection shapes.\n"); + } + } + +// return a const (!!) reference to the upwind + return *const_cast*>(m_spConvShape.get()); +} + +//////////////////////////////////////////////////////////////////////////////// +// register assemble functions +//////////////////////////////////////////////////////////////////////////////// + +#ifdef UG_DIM_1 +template<> +void ConvectionDiffusionFV1_cutElem:: +register_all_funcs(bool bHang) +{ + register_func > >(); + +/* +// switch assemble functions + if(!bHang) + { + register_func >(); + } + else + { + register_func >(); + } + */ +} +#endif + +#ifdef UG_DIM_2 +template<> +void ConvectionDiffusionFV1_cutElem:: +register_all_funcs(bool bHang) +{ + register_func > >(); + +/* +// switch assemble functions + if(!bHang) + { + register_func >(); + register_func >(); + register_func >(); + } + else + { + register_func >(); + register_func >(); + register_func >(); + } + */ +} +#endif + +#ifdef UG_DIM_3 +template<> +void ConvectionDiffusionFV1_cutElem:: +register_all_funcs(bool bHang) +{ + register_func > >(); + + /* + +// switch assemble functions + if(!bHang) + { + register_func >(); + register_func >(); + register_func >(); + register_func >(); + register_func >(); + register_func >(); + register_func >(); + register_func >(); + } + else + { + register_func >(); + register_func >(); + register_func >(); + register_func >(); + register_func >(); + register_func >(); + register_func >(); + register_func >(); + } +*/ +} +#endif + +template +template +void ConvectionDiffusionFV1_cutElem:: +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_A_expl_elem_fct(id, &T::template add_def_A_expl_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 ConvectionDiffusionFV1_cutElem; +#endif +#ifdef UG_DIM_2 +template class ConvectionDiffusionFV1_cutElem; +#endif +#ifdef UG_DIM_3 +template class ConvectionDiffusionFV1_cutElem; +#endif + +} // end namespace ConvectionDiffusionPlugin +} // namespace ug + +#endif /* __H__UG__LIB_DISC__CONVECTION_DIFFUSION__CONVECTION_DIFFUSION_FV1_CUTELEM_IMPL */ \ No newline at end of file diff --git a/fv1_cutElem/convection_diffusion_fv1_cutElem.h b/fv1_cutElem/convection_diffusion_fv1_cutElem.h new file mode 100644 index 0000000..0e77061 --- /dev/null +++ b/fv1_cutElem/convection_diffusion_fv1_cutElem.h @@ -0,0 +1,289 @@ +/* + * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt + * Author: 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__LIB_DISC__CONVECTION_DIFFUSION__CONVECTION_DIFFUSION_FV1_CUTELEM_ +#define __H__UG__LIB_DISC__CONVECTION_DIFFUSION__CONVECTION_DIFFUSION_FV1_CUTELEM_ + + +// ug4 headers +#include "lib_disc/spatial_disc/disc_util/conv_shape_interface.h" + +// plugin's internal headers +#include "../convection_diffusion_base.h" + + +#include "lib_disc/spatial_disc/disc_util/fv1Cut_geom.h" +#include "lib_disc/spatial_disc/disc_util/geom_provider.h" + +namespace ug{ +namespace ConvectionDiffusionPlugin{ + +// \ingroup lib_disc_elem_disc +/// \addtogroup convection_diffusion +/// \{ + +/// Discretization for the Convection-Diffusion Equation +/** + * This class implements the IElemDisc interface to provide element local + * assemblings for the convection diffusion equation. + * The Equation has the form + * \f[ + * \partial_t (m1*c + m2) - \nabla \left( D \nabla c - \vec{v} c - \vec{F} \right ) + * + r1 \cdot c + r2 = f + f2 + * \f] + * with + *
    + *
  • \f$ c \f$ is the unknown solution + *
  • \f$ m1 \equiv m(\vec{x},t) \f$ is the Mass Scaling Term + *
  • \f$ m2 \equiv m(\vec{x},t) \f$ is the Mass Term + *
  • \f$ D \equiv D(\vec{x},t) \f$ is the Diffusion Tensor + *
  • \f$ v \equiv \vec{v}(\vec{x},t) \f$ is the Velocity Field + *
  • \f$ F \equiv \vec{F}(\vec{x},t) \f$ is the Flux + *
  • \f$ r1 \equiv r(\vec{x},t) \f$ is the Reaction Rate + *
  • \f$ r2 \equiv r(\vec{x},t) \f$ is a Reaction Term + *
  • \f$ f \equiv f(\vec{x},t) \f$ is a Source Term + *
  • \f$ f2 \equiv f_2(\vec{x},t) \f$ is a Vector Source Term + *
+ * + * \tparam TDomain Domain + */ +template< typename TDomain> +class ConvectionDiffusionFV1_cutElem : public ConvectionDiffusionBase +{ + private: + /// Base class type + typedef ConvectionDiffusionBase base_type; + + /// Own type + typedef ConvectionDiffusionFV1_cutElem this_type; + + /// error estimator type + typedef SideAndElemErrEstData err_est_type; + + public: + /// World dimension + static const int dim = base_type::dim; + + public: + /// Constructor + ConvectionDiffusionFV1_cutElem(const char* functions, const char* subsets); + + /// set the upwind method + /** + * This method sets the upwind method used to upwind the convection. + * + * \param shapes upwind method + */ + void set_upwind(SmartPtr > shapes); + + private: + /// prepares assembling + virtual void prep_assemble_loop(); + + /// prepares the loop over all elements + /** + * This method prepares the loop over all elements. It resizes the Position + * array for the corner coordinates and schedules the local ip positions + * at the data imports. + */ + template + void prep_elem_loop(const ReferenceObjectID roid, const int si); + + /// prepares the element for assembling + /** + * This methods prepares an element for the assembling. The Positions of + * the Element Corners are read and the Finite Volume Geometry is updated. + * The global ip positions are scheduled at the data imports. + */ + template + void prep_elem(const LocalVector& u, GridObject* elem, const ReferenceObjectID roid, const MathVector vCornerCoords[]); + + /// finishes the loop over all elements + template + void fsh_elem_loop(); + + /// wrapper method for assembly of the local stiffness matrix using a finite volume scheme + template + void add_jac_A_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]); + + /// adapts local data for cut element assemling and calls 'add_jac_A_elem_local_local()' and 'add_jac_A_elem_local_boundary()' + /// for the assembly of the local stiffness matrix + template + void add_jac_A_elem_local_cut(TFVGeom& geo, const LocalVector& u, GridObject* elem, + const MathVector vCornerCoords[]); + + /// assembles the local stiffness matrix for locally adapted data due to a cut element + template + void add_jac_A_elem_local(TFVGeom& geo, const LocalVector& u, LocalMatrix& J, LocalVector& locU, GridObject* elem, + const MathVector vCornerCoords[]); + + /// assembles the boundary condition on a cut element + template + void add_jac_A_elem_boundary(TFVGeom& geo, LocalMatrix& J, const LocalVector& u, GridObject* elem, + const MathVector vCornerCoords[]); + + + /// assembles the local mass matrix using a finite volume scheme + template + void add_jac_M_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]); + + + /// wrapper method for assembly of the local defect using a finite volume scheme + template + void add_def_A_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]); + + /// adapts local data for cut element assembling, calls 'add_def_A_elem_local_local()', 'add_def_A_elem_local_boundary()' + /// and 'add_l2error_A_elem()' for the assembly of the local defect + template + void add_def_A_elem_local_cut(TFVGeom& geo, const LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]); + + /// assembles the local stiffness matrix for locally adapted data due to a cut element + template + void add_def_A_elem_local(TFVGeom& geo, const LocalVector& u, LocalVector& d, LocalVector& locU, GridObject* elem, + const MathVector vCornerCoords[]); + /// assembles the boundary condition on a cut element + template + void add_def_A_elem_boundary(TFVGeom& geo, LocalVector& d, LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]); + + /// computes the l2 error on each element => also on the 2 parts of a cut element + /// --> used for an interface-adapted l2 error computation + template + number add_l2error_A_elem(TFVGeom& geo, ReferenceObjectID roid, LocalVector& d, const LocalVector& u, GridObject* elem); + + + /// assembles the stiffness part of the local defect explicit reaction, reaction_rate and source + template + void add_def_A_expl_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]); + + /// assembles the mass part of the local defect + template + void add_def_M_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector vCornerCoords[]); + + /// assembles the local right hand side + template + void add_rhs_elem(LocalVector& d, GridObject* elem, const MathVector vCornerCoords[]); + template + void add_rhs_elem_local(TFVGeom& geo, LocalVector& d, GridObject* elem, const MathVector vCornerCoords[]); + + /// helper function to prepare local data due to a cut element for 'add_def_A_elem_local()' and 'add_jac_A_elem_local()' + template + void get_local_data(TFVGeom& geo, const LocalVector& u, const LocalIndices& ind, LocalVector& locUOut, MathMatrix& diffusionOut, LocalVector& jumpOut, LocalVector& jump_gradOut, LocalVector& sourceOut); + + public: + /// flag to set the analytical solution for the l2error computation within 'add_l2error_A_elem()' + size_t get_testCase() { return m_testCase; } + void set_testCase(size_t testCase) { m_testCase = testCase; } + + private: + /// abbreviation for the local solution + static const size_t _C_ = 0; + + using base_type::m_imDiffusion; + using base_type::m_imVelocity; + using base_type::m_imFlux; + using base_type::m_imSource; + using base_type::m_imSourceExpl; + using base_type::m_imVectorSource; + using base_type::m_imReactionRate; + using base_type::m_imReactionRateExpl; + using base_type::m_imReaction; + using base_type::m_imReactionExpl; + using base_type::m_imMassScale; + using base_type::m_imMass; + + using base_type::m_exGrad; + using base_type::m_exValue; + + protected: + /// method to compute the upwind shapes + SmartPtr > m_spConvShape; + + /// returns the updated convection shapes + typedef IConvectionShapes conv_shape_type; + const IConvectionShapes& get_updated_conv_shapes(const FVGeometryBase& geo); + + public: + /// type of trial space for each function used + virtual void prepare_setting(const std::vector& vLfeID, bool bNonRegularGrid); + + /// returns if hanging nodes are needed + virtual bool use_hanging() const; + + protected: + /// current regular grid flag + bool m_bNonRegularGrid; + + /// current shape function set (needed for GeomProvider::get()) + LFEID m_LFEID; + + /// flag to set the analytical solution for the l2error computation within 'add_l2error_A_elem()' + size_t m_testCase; + + /// register utils + /// \{ + void register_all_funcs(bool bHang); + template void register_func(); + /// \} + + private: + /// struct holding values of shape functions in IPs + struct ShapeValues + { + public: + void resize(std::size_t nEip, std::size_t nSip, std::size_t _nSh) + { + nSh = _nSh; + elemVals.resize(nEip); + sideVals.resize(nSip); + for (std::size_t i = 0; i < nEip; i++) elemVals[i].resize(nSh); + for (std::size_t i = 0; i < nSip; i++) sideVals[i].resize(nSh); + } + number& shapeAtElemIP(std::size_t sh, std::size_t ip) {return elemVals[ip][sh];} + number& shapeAtSideIP(std::size_t sh, std::size_t ip) {return sideVals[ip][sh];} + number* shapesAtElemIP(std::size_t ip) {return &elemVals[ip][0];} + number* shapesAtSideIP(std::size_t ip) {return &sideVals[ip][0];} + std::size_t num_sh() {return nSh;} + private: + std::size_t nSh; + std::vector > elemVals; + std::vector > sideVals; + } m_shapeValues; +}; + +// end group convection_diffusion +/// \} + +} // end ConvectionDiffusionPlugin +} // end namespace ug + + +#endif /*__H__UG__LIB_DISC__CONVECTION_DIFFUSION__CONVECTION_DIFFUSION_FV1_CUTELEM_*/ diff --git a/fv1_cutElem/diffusion_interface/diffusion_interface.h b/fv1_cutElem/diffusion_interface/diffusion_interface.h new file mode 100644 index 0000000..0eda600 --- /dev/null +++ b/fv1_cutElem/diffusion_interface/diffusion_interface.h @@ -0,0 +1,181 @@ +/* + * diffusion_interface.h + * + * Created on: 24.08.2017 + * Author: suze + */ + +#ifndef IMMERSED_INTERFACE_DIFFUSION_H_ +#define IMMERSED_INTERFACE_DIFFUSION_H_ + + +#ifdef UG_PARALLEL + #include "lib_grid/parallelization/load_balancer_util.h" +#endif + +#include "../convection_diffusion_fv1_cutElem.h" +#include "../../convection_diffusion_base.h" +#include "lib_disc/spatial_disc/immersed_util/interface_handler/interface_handler_two_sided_cut/interface_handler_diffusion.h" +#include "loc_to_glob_mapper_diffusion.h" + +namespace ug{ +namespace ConvectionDiffusionPlugin{ + + + +template < typename TDomain, typename TAlgebra> +class ImmersedInterfaceDiffusion + : public IImmersedInterface +{ + public: + /// world Dimension + static const int dim = TDomain::dim; + + /// Algebra type + typedef TAlgebra algebra_type; + + /// Type of algebra matrix + typedef typename algebra_type::matrix_type matrix_type; + + /// Type of algebra vector + typedef typename algebra_type::vector_type vector_type; + + typedef typename domain_traits::grid_base_object grid_base_object; + + ImmersedInterfaceDiffusion( + SmartPtr > ass, + SmartPtr > spMaster, + SmartPtr > interfaceProvider, + SmartPtr > cutElementHandler); + + // destructor + virtual ~ImmersedInterfaceDiffusion(){}; + + ////////////////////////////////////////////////////////////////////////////////////////// + // main methods + ////////////////////////////////////////////////////////////////////////////////////////// + + // general initialisation of set up data; + // most important: call of + // (1) 'initialize_interface()' and + // (2) 'update_interface_data()' to mark the cut elements and interface vertices + void init(vector_type& u, SmartPtr > spApproxSpace, const int baseLevel, + const int topLevel, bool bScaleDoFs); + + // mainly updates the BoolMarker for elements and vertices + void update(vector_type& u, SmartPtr > spApproxSpace, + const int baseLevel, const int topLevel, const number time); + + ////////////////////////////////////////////////////////////////////////////////////////// + // helper methods for init() and update() + ////////////////////////////////////////////////////////////////////////////////////////// + + ////////////////////////////////////////////////////////////////////////////////////////// + /// Info - 'initialize_interface()': Main infos: see CutElementHandler + /// Important for call here: + /// --> counts number of cut elements! + /// --> called during init() + ////////////////////////////////////////////////////////////////////////////////////////// + + const size_t initialize_interface(ConstSmartPtr dd) + { return m_spCutElementHandler->initialize_interface(dd); } + const size_t initialize_interface_Nitsche(vector_type& u, ConstSmartPtr dd); + + ////////////////////////////////////////////////////////////////////////////////////////// + /// methods for writing the source term and boundary conditions on the interface + ////////////////////////////////////////////////////////////////////////////////////////// + + void set_source_data_lua(const number interfaceSourceVal) + { m_spInterfaceHandlerLocal->set_source_data_lua(interfaceSourceVal); } + void set_jump_data_lua(const number interfaceJumpVal) + { m_spInterfaceHandlerLocal->set_jump_data_lua(interfaceJumpVal); } + void set_jump_grad_data_lua(const MathVector<2>& interfaceJumpGradVal) + { m_spInterfaceHandlerLocal->set_jump_grad_data_lua(interfaceJumpGradVal); } + void set_diffusion_data_lua(const MathVector<2>& diffusionCoeffs) + { m_spInterfaceHandlerLocal->set_diffusion_coeff_data_lua(diffusionCoeffs); } + + ////////////////////////////////////////////////////////////////////////////////////////// + /// lua-methods for set up: + ////////////////////////////////////////////////////////////////////////////////////////// + + // the 'threshold' defines the bandwidth around the immersed interface, in which a node + // counts as 'OUTSIDE' or 'ON_INTERFACE' during call of 'CutElementHandler::is_outside()', + // 'CutElementHandler::is_nearInterface() + void initialize_threshold(TDomain& domain, const int baseLevel, const int topLevel); + void set_threshold(size_t level, const number threshold) + { m_spCutElementHandler->set_threshold(level, threshold); } + number MeanElementDiameter(TDomain& domain, int level); + + // If StdFV-assembling is ON, NO new 'vCornerCoords' will be computed on the cut elements. + // The original nodes ACCROSS the interface (and ON the euclidian mesh) will be chosen for the + // computation of the solution of the interface + // ==> the standard shape functions as in common 'ficticious domain' methods will be used + // ==> the shape functions will NOT be 1 ON the interface and the gradient will NOT + // point normal to the interface + void set_StdFV_assembling(bool bValue) { m_spInterfaceHandlerLocal->set_StdFV_assembling(bValue);} + bool StdFV_assembling() { return m_spInterfaceHandlerLocal->StdFV_assembling(); } + + + ////////////////////////////////////////////////////////////////////////////////////////// + /// lua-methods for output + ////////////////////////////////////////////////////////////////////////////////////////// + + void set_analytic_solution(vector_type& u, SmartPtr > spApproxSpace, + SmartPtr mg, const int topLevel); + + // adjustment of solution vector in order to compute the error WITHOUT nodes near the interface: + // ==> (1) remove additional nodes ON the interface, not lying on the original grid + // (2) set solution in nodes lying NEAR the interface, but IN the original grid to zero + void adjust_for_error(vector_type& u, vector_type& uCopy, SmartPtr > spApproxSpace, + SmartPtr mg, const int topLevel); + + double compute_solution_value(const MathVector& vrtPos); + + // returns the interface adapted l2 error, computed during 'ConvectionDiffusionFV1_cutElem::add_def_A_elem()': + number get_L2Error() { return m_spInterfaceHandlerLocal->get_L2Error(); } + + // returns the number of DoFs on the original grid + size_t get_numDoFs() { return m_spInterfaceMapper->get_numDoFs(); } + + // boolian used for perform the Nitsche-method (i.e. CutElem-method) for the treatment of the immersed interface + void set_Nitsche(bool bNitsche) { return m_spInterfaceHandlerLocal->set_Nitsche(bNitsche); } + + // setting and getting flag for printing of cut-element data into file: + void set_print_cutElemData(bool bValue) { m_spInterfaceHandlerLocal->set_print_cutElemData(bValue); } + + // returns the number of DoFs on the original grid + size_t get_numCutElements(const int gridlevel, const size_t prtIndex) + { + ConstSmartPtr dd = m_spApproxSpace->dof_distribution(GridLevel(gridlevel, GridLevel::LEVEL)); + const int levIndex = m_spCutElementHandler->get_Index(gridlevel, dd); + + return m_spCutElementHandler->get_numCutElements(levIndex); + } + + ////////////////////////////////////////////////////////////////////////////////////////// + /// class member + ////////////////////////////////////////////////////////////////////////////////////////// + + private: + /// current ApproxSpace + SmartPtr > m_spApproxSpace; + + SmartPtr > m_spInterfaceProvider; + SmartPtr > m_spCutElementHandler; + SmartPtr > m_spInterfaceHandlerLocal; + + SmartPtr > m_spInterfaceMapper; + +}; + +} // end namespace ConvectionDiffusionPlugin +} // end namespace ug + + +#include "diffusion_interface_impl.h" + +#endif /* IMMERSED_INTERFACE_DIFFUSION_H_ */ + + + + diff --git a/fv1_cutElem/diffusion_interface/diffusion_interface_impl.h b/fv1_cutElem/diffusion_interface/diffusion_interface_impl.h new file mode 100644 index 0000000..baf1895 --- /dev/null +++ b/fv1_cutElem/diffusion_interface/diffusion_interface_impl.h @@ -0,0 +1,355 @@ + +/* + * diffusion_interface_impl.h + * + * Created on: 24.08.2017 + * Author: suze + */ + +#ifndef IMMERSED_INTERFACE_DIFFUSION_IMPL_H_ +#define IMMERSED_INTERFACE_DIFFUSION_IMPL_H_ + +#include "lib_disc/spatial_disc/disc_util/geom_provider.h" + +namespace ug { +namespace ConvectionDiffusionPlugin { + +/////////////////////////////////////////////////////////// +// Implementation of the methods class +// 'ImmersedInterfaceDiffusion' +/////////////////////////////////////////////////////////// + + +template +ImmersedInterfaceDiffusion::ImmersedInterfaceDiffusion( + SmartPtr > ass, + SmartPtr > spMaster, + SmartPtr > interfaceProvider, + SmartPtr > cutElementHandler) : + m_spInterfaceProvider(interfaceProvider), + m_spCutElementHandler(cutElementHandler), + m_spInterfaceHandlerLocal(new InterfaceHandlerLocalDiffusion(interfaceProvider, cutElementHandler)), + m_spInterfaceMapper(new DiffusionInterfaceMapper (m_spInterfaceHandlerLocal)) +{ + if (interfaceProvider->num_particles() == 0) + UG_THROW("ImmersedInterfaceDiffusion::Constructor(): no particles initializen in 'globalHandler\n"); + +// initialize singleton and set local handler + typedef DimFV1CutGeometry > TFVGeom; + TFVGeom& geo = GeomProvider::get(LFEID(LFEID::LAGRANGE, dim, 1),1); + geo.set_interface_handler(m_spInterfaceHandlerLocal); + +// initialize mapper within domainDisc: + SmartPtr > assAdapt = ass->ass_tuner(); + assAdapt->set_mapping(m_spInterfaceMapper.get()); + +// needs to be enabled, in order to call 'spAssTuner->modify_LocalData()' during element disc assembling + assAdapt->enable_modify_solution(true); + + +} + + +template +void ImmersedInterfaceDiffusion:: +init(vector_type& u, SmartPtr > spApproxSpace, const int baseLevel, + const int topLevel, bool bScaleDoFs) +{ + // get data + m_spApproxSpace = spApproxSpace; + ConstSmartPtr dd = spApproxSpace->dof_distribution(GridLevel(topLevel, GridLevel::LEVEL)); + + // get the level Index ONLY for the toplevel in order to resize the data for the toplevel + // --> not implemented for multigrid already! + const int levIndex = m_spCutElementHandler->get_Index(topLevel, dd); + + // initialize the orientation of the interface: = 1, i.e. inside the circle line is outside the domain + m_spCutElementHandler->set_orientation(1); + +// call of 'update_interface_data()' (via init() of CutElementHandler): +// sets the marker (INSIDE, OUTSIDE, CUT_BY_INTERFACE) to each element +// and fills the lists of cut elements + m_spCutElementHandler->template init(dd, baseLevel, topLevel); + +// store the number of regular DoFs (u.size()) and additional DoFs on the immersed interface: + size_t numDoFs = u.size(); + size_t num_newDoFs = m_spCutElementHandler->get_numCutElements(levIndex); + m_spInterfaceMapper->set_numDoFs(numDoFs); + m_spInterfaceMapper->set_numNewDoFs(num_newDoFs); + + + if ( m_spInterfaceHandlerLocal->get_Nitsche() ) + { + initialize_interface_Nitsche(u, dd); + const size_t num_NitscheDoFs = m_spInterfaceHandlerLocal->get_num_NitscheDoFs(); + + num_newDoFs = num_NitscheDoFs; + m_spInterfaceMapper->set_numNewDoFs(num_newDoFs); + + } + +// initialize the scaling of new DoFs: *2 for jumping values (strong discontinuity) at interface +// => 1 DoFs is defined for the solution on each side of the interface +// values for new DoFs are set to 0.0 by the 'resize()'-method (see vector.h): + if ( bScaleDoFs ) u.resize(numDoFs + 2*num_newDoFs); + else u.resize(numDoFs + num_newDoFs); + + m_spInterfaceMapper->set_bScaleDoFs(bScaleDoFs); + m_spInterfaceHandlerLocal->set_bScaleDoFs(bScaleDoFs); + +// initialize the value for the l2 error computaton + m_spInterfaceHandlerLocal->L2Error_init(); + +} + +template +void ImmersedInterfaceDiffusion:: +update(vector_type& u, SmartPtr > spApproxSpace, + const int baseLevel, const int topLevel, const number time) +{ +// get data + ConstSmartPtr dd = spApproxSpace->dof_distribution(GridLevel(topLevel, GridLevel::LEVEL)); + int topLev = spApproxSpace->num_levels()-1; + if ( topLev != topLevel ) + UG_THROW("update: parameter 'topLevel' = " << topLevel << " != " + << topLev << "current top leven! \n"); + +// update data: update BoolMarker and cut element lists + m_spCutElementHandler->template init(dd, baseLevel, topLevel); + +} + +template +double ImmersedInterfaceDiffusion:: +compute_solution_value(const MathVector& vrtPos) +{ + double kappa_1 = 1.0; + double kappa_2 = 10.0; + double sqR = 0.4*0.4; + double dist_x = vrtPos[0] - 0.1; + double dist_y = vrtPos[1] - 0.2; + double sqDist = dist_x*dist_x+dist_y*dist_y; + + double value = -4*kappa_1*kappa_2*kappa_2*sqR*sqDist + 2*sqR*sqR*kappa_2*(2*kappa_1*kappa_2 - 1); + + double dist = sqrt(sqDist); + if ( dist > 0.4 ) + value = -2*kappa_2*sqDist*sqDist; + + return value; +} + +template +void ImmersedInterfaceDiffusion:: +set_analytic_solution(vector_type& u, SmartPtr > spApproxSpace, SmartPtr mg, const int topLevel) +{ + ConstSmartPtr dd = spApproxSpace->dof_distribution(GridLevel(topLevel, GridLevel::LEVEL)); + + typedef MathVector position_type; + typedef Attachment position_attachment_type; + typedef Grid::VertexAttachmentAccessor position_accessor_type; + + position_attachment_type m_aPos = GetDefaultPositionAttachment(); + position_accessor_type m_aaPos; + + if(!mg->has_attachment(m_aPos)) + mg->attach_to(m_aPos); + m_aaPos.access(*mg, m_aPos); + + typedef typename domain_traits::grid_base_object grid_base_object; + + typename DoFDistribution::traits::const_iterator iter, iterEnd; + iter = dd->template begin(); + iterEnd = dd->template end(); + + // loop elements in order to compute the volume and set rhs: + for( ; iter != iterEnd; iter++) + { + // get element + grid_base_object* elem = *iter; + std::vector > vCornerCoords; + CollectCornerCoordinates(vCornerCoords, *elem, m_aaPos); + + std::vector ind; + dd->dof_indices(elem, 0, ind); + + // loop vertices + for(size_t i = 0; i < elem->num_vertices(); ++i) + { + double value = compute_solution_value(vCornerCoords[i]); + DoFRef(u, ind[i]) = value; + } + } + +} + +template +void ImmersedInterfaceDiffusion:: +adjust_for_error(vector_type& u, vector_type& uCopy, SmartPtr > spApproxSpace, + SmartPtr mg, const int topLevel) +{ +// (1) remove additional nodes ON the interface, not lying on the original grid +// --> resize vector 'u' to original size, i.e. decrease size + size_t numDoFsCopy = uCopy.size(); + + u.resize(numDoFsCopy); + + +// get data + typedef MathVector position_type; + typedef Attachment position_attachment_type; + typedef Grid::VertexAttachmentAccessor position_accessor_type; + + position_attachment_type m_aPos = GetDefaultPositionAttachment(); + position_accessor_type m_aaPos; + + if(!mg->has_attachment(m_aPos)) mg->attach_to(m_aPos); + m_aaPos.access(*mg, m_aPos); + + ConstSmartPtr dd = spApproxSpace->dof_distribution(GridLevel(topLevel, GridLevel::LEVEL)); + + typedef typename domain_traits::grid_base_object grid_base_object; + typename DoFDistribution::traits::const_iterator iter, iterEnd; + iter = dd->template begin(); + iterEnd = dd->template end(); + +// (2) set solution to zero in near-interface nodes +// --> loop elements in order to set solution to zero on cut element nodes of vector 'u': + for( ; iter != iterEnd; iter++) + { + // get element + grid_base_object* elem = *iter; + + ElementModus elemModus = m_spInterfaceHandlerLocal->compute_element_modus(elem); + + // choose only cut elements for adjustment + bool do_adjust = false; + switch(elemModus) + { + case INSIDE_DOM: break; + case OUTSIDE_DOM: break; + + case CUT_BY_INTERFACE: do_adjust = true; + break; + default: + throw(UGError("Error in InterfaceHandlerLocalDiffusion::update(): switch(m_elemModus)!")); + } + + if ( !do_adjust ) continue; + + // get local indices + std::vector ind; + dd->dof_indices(elem, 0, ind); + + // loop vertices and set solution to zero + for(size_t i = 0; i < elem->num_vertices(); ++i) + DoFRef(u, ind[i]) = 0.0; + + } + +} + +template +const size_t ImmersedInterfaceDiffusion:: +initialize_interface_Nitsche(vector_type& u, ConstSmartPtr dd) +{ + m_spInterfaceHandlerLocal->m_MapInserted_Nitsche.clear(); + + typedef typename domain_traits::grid_base_object grid_base_object; + + typename DoFDistribution::traits::const_iterator iter, iterEnd; + iter = dd->template begin(); + iterEnd = dd->template end(); + + size_t num_cutElements = 0; + + // loop elements in order to compute the volume and set rhs: + for( ; iter != iterEnd; iter++) + { + // get element + grid_base_object* elem = *iter; + + ElementModus elemModus = m_spCutElementHandler->compute_element_modus(elem); + + if ( elemModus == CUT_BY_INTERFACE ) + { + num_cutElements += 1; + + // get global indices + LocalIndices ind; + dd->indices(elem, ind, false); + + // fill 'm_MapInserted_Nitsche' with global indices: + for(size_t i = 0; i < 3; ++i) + m_spInterfaceHandlerLocal->get_or_insert_indexPair_Nitsche(ind.index(0, i)); + + } + } + + return num_cutElements; + +} + + +template +number ImmersedInterfaceDiffusion:: +MeanElementDiameter(TDomain& domain, int level) +{ + typedef typename domain_traits::grid_base_object TElem; + typedef typename geometry_traits::iterator ListIter; + + ListIter iter = domain.grid()->template begin(level); + ListIter iterEnd = domain.grid()->template end(level); + + number mean = 0.0; + size_t numIter = 0; + for (; iter != iterEnd; ++iter) { + mean += ElementDiameterSq(*domain.grid(), domain.position_accessor(), + *iter); + numIter++; + } + + mean = mean / numIter; + +#ifdef UG_PARALLEL +// share value between all procs + pcl::ProcessCommunicator com; +// ToDO: PCL_RO_MIN oder MAX oder doch mean noch berechnen m.H.v. /numProcs??? +//UG_THROW("in MeanElementDiameter(): was macht 'allredue' im Parallelen???\n"); + mean = com.allreduce(mean, PCL_RO_MIN); +#endif + + return std::sqrt(mean); +} + +template +void ImmersedInterfaceDiffusion:: +initialize_threshold(TDomain& domain, const int baseLevel, const int topLevel) +{ + UG_THROW("initialize_threshold(): not tested for diffusion case! \n"); + + if (baseLevel < 0) + UG_THROW( + "initialize_threshold(): no cast of baselevel from 'int' tp 'size_t' possible! \n"); + if (topLevel < 0) + UG_THROW( + "initialize_threshold(): no cast of toplevel from 'int' tp 'size_t' possible! \n"); + +// compute level-dependent value for threshold: + for (size_t lev = baseLevel; lev <= (size_t) topLevel; ++lev) + { + //const number maxLength = MaxElementDiameter(domain, lev); + const number meanLength = MeanElementDiameter(domain, lev); + +// set_threshold(lev, meanLength * meanLength); + set_threshold(lev, 0.25 * meanLength*meanLength); + } + + +} + + +} // end namespace ConvectionDiffusionPlugin +} // end namespace ug + +#endif /* IMMERSED_INTERFACE_DIFFUSION_IMPL_H_ */ diff --git a/fv1_cutElem/diffusion_interface/loc_to_glob_mapper_diffusion.h b/fv1_cutElem/diffusion_interface/loc_to_glob_mapper_diffusion.h new file mode 100644 index 0000000..42c0f6c --- /dev/null +++ b/fv1_cutElem/diffusion_interface/loc_to_glob_mapper_diffusion.h @@ -0,0 +1,168 @@ +/* + * diffusion_interface_handler_local.h + * + * Created on: 19.01.2015 + * Author: suze + */ + +#ifndef DIFFUSION_INTERFACE_MAPPER_H_ +#define DIFFUSION_INTERFACE_MAPPER_H_ + +#include "lib_disc/spatial_disc/immersed_util/immersed_interface_base.h" + +namespace ug{ +namespace ConvectionDiffusionPlugin{ + + +template +class DiffusionInterfaceMapper : public IInterfaceMapper +{ + + public: + /// World dimension + static const int dim = TDomain::dim; + /// Algebra type + typedef TAlgebra algebra_type; + + /// Type of algebra matrix + typedef typename algebra_type::matrix_type matrix_type; + + /// Type of algebra vector + typedef typename algebra_type::vector_type vector_type; + + /// Type of geometric base object + typedef typename domain_traits::grid_base_object grid_base_object; + + DiffusionInterfaceMapper(){}; + + DiffusionInterfaceMapper(SmartPtr > localHandler) + : m_numGridNodes(0), + m_numNewDoFs(0), + m_resized(false), + m_resized_defect(false), + m_scaleDoFs(false), + m_spInterfaceHandlerLocal(localHandler) + {} + + virtual ~DiffusionInterfaceMapper() {} + + /////////////////////////////////////////////////////////////////////////////// + /// + /// base class methods not needed for 'DiffusionInterfaceMapper' class + /// + /////////////////////////////////////////////////////////////////////////////// + + /// modifies local solution vector for adapted defect computation + void modify_LocalData(LocalMatrix& locJ, LocalVector& locU, + ConstSmartPtr dd){}; + void modify_LocalData(LocalVectorTimeSeries& uT, LocalMatrix& locJ, LocalVector& locU, + ConstSmartPtr dd){}; + + void modify_LocalData(LocalVector& locD, LocalVector& tmpLocD, LocalVector& locU, + ConstSmartPtr dd){}; + void modify_LocalData(LocalVectorTimeSeries& uT, LocalVector& locD, LocalVector& tmpLocD, + LocalVector& locU, ConstSmartPtr dd, size_t t){}; + + /////////////////////////////////////////////////////////////////////////////// + /// + /// base class methods and helper methods called by them + /// + /////////////////////////////////////////////////////////////////////////////// + + /// base method: send local entries to global rhs + void add_local_vec_to_global(vector_type& vec, const LocalVector& lvec, + ConstSmartPtr dd); + /// methods called by base method for cut element case: special mapping due to new DoFs! + void add_local_vec_to_global_interface(vector_type& vec, const LocalVector& lvec, + ConstSmartPtr dd); + /// methods called by base method for Nitsche-treatment on cut elements + void add_local_vec_to_global_Nitsche(vector_type& vec, const LocalVector& lvec, + ConstSmartPtr dd); + + /// base method: send local entries to global matrix + void add_local_mat_to_global(matrix_type& mat, const LocalMatrix& lmat, + ConstSmartPtr dd); + /// methods called by base method for cut element case: special mapping due to new DoFs! + void add_local_mat_to_global_interface(matrix_type& mat, const LocalMatrix& lmat, + ConstSmartPtr dd); + void add_local_mat_to_global_Nitsche(matrix_type& mat, const LocalMatrix& lmat, + ConstSmartPtr dd); + + /// sets all non-DoFs to identity rows + /// --> or diffusion, iff only INSIDE circle-computation only! ==> all other DoFs set to Dirichlet! + void adjust_mat_global(matrix_type& mat, const LocalMatrix& lmat, ConstSmartPtr dd); + + + /////////////////////////////////////////////////////////////////////////////// + /// REMARK: + /// During DomainDiscretization::assemble_jacobian: + /// calling + /// ---> m_spAssTuner->modify_GlobalSol(pModifyMemory, vSol, dd); + /////////////////////////////////////////////////////////////////////////////// + + /// instead of modifying global data: the computed values at the interface DoFs get + /// written/stored into data of class 'InterfaceHandlerLocalDiffusion::m_verticesValue' + /// via call of 'set_interface_values()' (for each call of domainDisc, i.e. newton step) + + void modify_GlobalSol(vector_type& vecMod, const vector_type& vec, + ConstSmartPtr dd); + + void modify_GlobalSol(SmartPtr > vSolMod, + ConstSmartPtr > vSol, + ConstSmartPtr dd){}; + + /////////////////////////////////////////////////////////////// + /// new methods: + /////////////////////////////////////////////////////////////// + + /// sets dirichlet rows for non-DoFs (not yet needed here) + void set_identity_mat_constraint(matrix_type& mat, const LocalMatrix& lmat, + ConstSmartPtr dd); + + // access methods to local data stored in class 'InterfaceHandlerLocalDiffusion' + // REMARK: the contribution of the boundary condition on the immersed interface + // were assembled within the ElemDisc 'ConvectionDiffusionFV1_cutElem' + // and stored in 'InterfaceHandlerLocalDiffusion' to access it here for + // performing the local-to-global mapping + LocalMatrix& get_local_jacobian(const ReferenceObjectID roid) + { return m_spInterfaceHandlerLocal->get_local_jacobian(roid); } + LocalVector& get_local_defect(const ReferenceObjectID roid) + { return m_spInterfaceHandlerLocal->get_local_defect(roid); } + + + + // writes the solution in the interface nodes, i.e. the NEW DoFs, + // into data of 'InterfaceHandlerLocalDiffusion' + void set_interface_values(const std::vector verticesValues) + { m_spInterfaceHandlerLocal->set_interface_values(verticesValues); } + + // called during 'InterfaceHandlerLocalDiffusion::init() + void set_numDoFs(const size_t numDoFs) { m_numGridNodes = numDoFs;} + size_t get_numDoFs() { return m_numGridNodes;} + void set_numNewDoFs(const size_t numNewDoFs) { m_numNewDoFs = numNewDoFs;} + + // returns the value, by which the global matrix and vector needs to be resized + // REMARK: for a jump in the solution, 2 new DoFs will be placed for each interface node + const size_t get_resize_measure(const size_t numDoFs); + + void set_bScaleDoFs(bool bScaleDoF) { m_scaleDoFs = bScaleDoF; } + + private: + size_t m_numGridNodes; // number of grid nodes, i.e. of DoFs WITHOUT interface nodes + size_t m_numNewDoFs; // number of interface nodes, i.e. additional DoFs + bool m_resized; + bool m_resized_defect; + bool m_scaleDoFs; // if m_scaleDoFs = true: 2 new DoFs will be placed for + // each interface node (for jump in value) + + SmartPtr > m_spInterfaceHandlerLocal; + +}; + +}// namespace ConvectionDiffusionPlugin +} // end ug namespace + +#include "loc_to_glob_mapper_diffusion_impl.h" + + +#endif /* DIFFUSION_INTERFACE_MAPPER_H_ */ diff --git a/fv1_cutElem/diffusion_interface/loc_to_glob_mapper_diffusion_impl.h b/fv1_cutElem/diffusion_interface/loc_to_glob_mapper_diffusion_impl.h new file mode 100644 index 0000000..ae0fd5e --- /dev/null +++ b/fv1_cutElem/diffusion_interface/loc_to_glob_mapper_diffusion_impl.h @@ -0,0 +1,517 @@ +/* + * diffusion_interface_handler_local_impl.h + * + * Created on: 15.01.2015 + * Author: suze + */ + +#ifndef DIFFUSION_INTERFACE_MAPPER_IMPL_H_ +#define DIFFUSION_INTERFACE_MAPPER_IMPL_H_ + + +namespace ug{ +namespace ConvectionDiffusionPlugin{ + + +template +void DiffusionInterfaceMapper:: +set_identity_mat_constraint(matrix_type& mat, const LocalMatrix& lmat, ConstSmartPtr dd) +{ + + const LocalIndices& rowInd = lmat.get_row_indices(); + + for (size_t fct1 = 0; fct1 < lmat.num_all_row_fct(); ++fct1) + for (size_t dof1 = 0; dof1 < lmat.num_all_row_dof(fct1); ++dof1) + { + const size_t rowIndex = rowInd.index(fct1, dof1); + const size_t rowComp = rowInd.comp(fct1, dof1); + + if ( fabs(lmat.value(fct1, dof1, fct1, dof1)) < 0.000001 + && BlockRef(mat(rowIndex, rowIndex), rowComp, rowComp) == 0.0 ) + BlockRef(mat(rowIndex, rowIndex), rowComp, rowComp) = 1.0; + } + +} + + +template +void DiffusionInterfaceMapper:: +adjust_mat_global(matrix_type& mat, const LocalMatrix& lmat, ConstSmartPtr dd) +{ + DoFIndex index; + + for (size_t i = 0; i < m_numGridNodes; ++i) + { + // set all entries, which were not assembled by the ElemDisc to the identity: + if ( BlockRef(mat(i, i), 0, 0) == 0.0 ) + BlockRef(mat(i, i), 0, 0) = 1.0; + } +} + + +template +void DiffusionInterfaceMapper:: +modify_GlobalSol(vector_type& vecMod, const vector_type& vec, ConstSmartPtr dd) +{ + size_t numDoFs = vecMod.size(); + const size_t numNewDoFs = numDoFs - m_numGridNodes; + + UG_LOG("---------------modify_GlobalSol--------------\n"); + UG_LOG(" vecMod.size(): " << numDoFs << "\n"); + UG_LOG(" m_numGridNodes: " << m_numGridNodes << "\n"); + UG_LOG(" computed numNewDoFs: " << numNewDoFs << "\n"); + UG_LOG(" m_numNewDoFs: " << m_numNewDoFs << "\n"); + + DoFIndex index; + std::vector verticesValue; + verticesValue.clear(); + + // loop all interface nodes, i.e new DoFs, and store the values in + // member 'm_verticesValue' of class 'InterfaceHandlerLocalDiffusion': + for (size_t i = 0; i < numNewDoFs; ++i) + { + // get GLOBAL index of the interface nodes + index = DoFIndex(m_numGridNodes + i,0); + // write value to local data + verticesValue.push_back(DoFRef(vec, index)); + } + +// write solution values to member 'm_verticesValue' of class 'InterfaceHandlerLocalDiffusion': + set_interface_values(verticesValue); + +} + +template +const size_t DiffusionInterfaceMapper:: +get_resize_measure(const size_t numDoFs) +{ + size_t numAllDoFs = m_numGridNodes + m_numNewDoFs; + + // if m_scaleDoFs = true: 2 new DoFs will be placed for each interface node + if ( m_scaleDoFs ) + numAllDoFs = m_numGridNodes + 2*m_numNewDoFs; + + return numAllDoFs - numDoFs; + +} + +template +void DiffusionInterfaceMapper:: +add_local_mat_to_global(matrix_type& mat, const LocalMatrix& lmat, ConstSmartPtr dd) +{ +// reset print modus for cut element data to 'false' +// => only printed during assemling of defect, i.e. ONE loop over all element + m_spInterfaceHandlerLocal->set_print_cutElemData(false); + + size_t numAllDoFs = m_numGridNodes + m_numNewDoFs; + +// resize global matrix dynamically: + const size_t diffDoFs = get_resize_measure(mat.num_rows()); + +// resize global defect only ONCE: + if ( diffDoFs > 0 ) + { + mat.resize_and_keep_values(numAllDoFs, numAllDoFs); + } + else if ( diffDoFs == 0 ) + { // no resizing + } + else if ( diffDoFs < 0 ) + { UG_THROW("error in add_local_mat_to_global: diffDofs < 0\n"); } + + + ElementModus modus = m_spInterfaceHandlerLocal->elementModus(); + + switch(modus) + { + case OUTSIDE_DOM: // call usual local-to-global-mapping + AddLocalMatrixToGlobal(mat, lmat); + break; + case INSIDE_DOM: // call usual local-to-global-mapping + AddLocalMatrixToGlobal(mat, lmat); + break; + case CUT_BY_INTERFACE: // call adapted local-to-global-mapping + if ( m_spInterfaceHandlerLocal->get_Nitsche() ) + add_local_mat_to_global_Nitsche(mat, lmat, dd); + else + add_local_mat_to_global_interface(mat, lmat, dd); + break; + default: + throw(UGError("Error in IInterfaceMapper::add_local_mat_to_global()!")); + + } + +} + + +template +void DiffusionInterfaceMapper:: +add_local_vec_to_global(vector_type& vec, const LocalVector& lvec, ConstSmartPtr dd) +{ + size_t numAllDoFs = m_numGridNodes + m_numNewDoFs; + +// resize global vector dynamically: + const size_t diffDoFs = get_resize_measure(vec.size()); + +// resize global defect ONCE: + if ( diffDoFs > 0 ) + { + vec.resize(numAllDoFs, true); + vec.set(0.0); + } + else if ( diffDoFs == 0 ) + { // no resizing + } + else if ( diffDoFs < 0 ) + { UG_THROW("error in add_local_vec_to_global: diffDofs < 0\n"); } + + + ElementModus modus = m_spInterfaceHandlerLocal->elementModus(); + + switch(modus) + { + case OUTSIDE_DOM: // call usual local-to-global-mapping + AddLocalVector(vec, lvec); + break; + case INSIDE_DOM: // call usual local-to-global-mapping + AddLocalVector(vec, lvec); + break; + case CUT_BY_INTERFACE: // call adapted local-to-global-mapping + if ( m_spInterfaceHandlerLocal->get_Nitsche() ) + add_local_vec_to_global_Nitsche(vec, lvec, dd); + else + add_local_vec_to_global_interface(vec, lvec, dd); + break; + default: + throw(UGError("Error in IInterfaceMapper::add_local_vec_to_global()!")); + + } +} + +template +void DiffusionInterfaceMapper:: +add_local_mat_to_global_Nitsche(matrix_type& mat, const LocalMatrix& lmat, ConstSmartPtr dd) +{ + const LocalIndices& rowInd = lmat.get_row_indices(); + const LocalIndices& colInd = lmat.get_col_indices(); + + DoFIndex indexRow, indexCol; + +/////////////////////////////////////////////////////////////// +/// FIRST: add loc to glob for locJ_tri: +/////////////////////////////////////////////////////////////// + + LocalMatrix locJ_tri = get_local_jacobian(ROID_TRIANGLE); + const bool shift_global_index_tri = m_spInterfaceHandlerLocal->get_index_shift_tri(); + + if ( lmat.num_all_row_dof(0) != locJ_tri.num_all_row_dof(0) ) + UG_THROW("in 'add_local_mat_to_global_Nitsche': non-consistent sizees!\n"); + + size_t numAllDoFs = m_numGridNodes; + + if ( shift_global_index_tri ) + numAllDoFs = m_numGridNodes + m_numNewDoFs; + + for (size_t dof1 = 0; dof1 < locJ_tri.num_all_row_dof(0); ++dof1) + { + size_t global_index1 = rowInd.index(0, dof1); + + // if global_index1 is index of m_vertex: compute global index: + if ( m_spInterfaceHandlerLocal->lies_onInterface(dof1) ) + global_index1 = numAllDoFs + m_spInterfaceHandlerLocal->get_index_for_global_index_Nitsche(global_index1); + + indexRow = DoFIndex(global_index1, rowInd.comp(0, dof1)); + + for (size_t dof2 = 0; dof2 < locJ_tri.num_all_col_dof(0); ++dof2) + { + size_t global_index2 = rowInd.index(0, dof2); + + // if global_index2 is index of m_vertex: compute global index: + if ( m_spInterfaceHandlerLocal->lies_onInterface(dof2) ) + global_index2 = numAllDoFs + m_spInterfaceHandlerLocal->get_index_for_global_index_Nitsche(global_index2); + + indexCol = DoFIndex(global_index2, colInd.comp(0, dof2)); + + // finally add loc to glob: + DoFRef(mat, indexRow, indexCol) += locJ_tri.value(0, dof1, 0, dof2); + + } + } + +/////////////////////////////////////////////////////////////// +/// SECOND: add loc to glob for locJ_tri: +/// -> same as FIRST, but: +/// (1) lmat instead of locJ_tri +/// (2) if ( !m_spInterfaceHandlerLocal->lies_onInterface(dof) ) +/// instead of +/// if ( m_spInterfaceHandlerLocal->lies_onInterface(dof) ) +/////////////////////////////////////////////////////////////// + + for (size_t dof1 = 0; dof1 < lmat.num_all_row_dof(0); ++dof1) + { + size_t global_index1 = rowInd.index(0, dof1); + + // if global_index1 is index of m_vertex: compute global index: + if ( !m_spInterfaceHandlerLocal->lies_onInterface(dof1) ) + global_index1 = numAllDoFs + m_spInterfaceHandlerLocal->get_index_for_global_index_Nitsche(global_index1); + + indexRow = DoFIndex(global_index1, rowInd.comp(0, dof1)); + + for (size_t dof2 = 0; dof2 < lmat.num_all_col_dof(0); ++dof2) + { + size_t global_index2 = rowInd.index(0, dof2); + + // if global_index2 is index of m_vertex: compute global index: + if ( !m_spInterfaceHandlerLocal->lies_onInterface(dof2) ) + global_index2 = numAllDoFs + m_spInterfaceHandlerLocal->get_index_for_global_index_Nitsche(global_index2); + + indexCol = DoFIndex(global_index2, colInd.comp(0, dof2)); + + // finally add loc to glob: + DoFRef(mat, indexRow, indexCol) += lmat.value(0, dof1, 0, dof2); + + } + } + + +} + +//get_real_index(dof): +// if dof NOT on interface: returns orig corner index +// if dof ON interface: returns index of m_vertex-array + +template +void DiffusionInterfaceMapper:: +add_local_mat_to_global_interface(matrix_type& mat, const LocalMatrix& lmat, ConstSmartPtr dd) +{ + const LocalIndices& rowInd = lmat.get_row_indices(); + const LocalIndices& colInd = lmat.get_col_indices(); + + DoFIndex indexRow, indexCol; + +///////////////////////////////////////////////////////////////////////////////// +/// FIRST: add the boundary contribution stored in 'locJ_tri' to global matrix: +///////////////////////////////////////////////////////////////////////////////// + + LocalMatrix locJ_tri = get_local_jacobian(ROID_TRIANGLE); + const bool shift_global_index_tri = m_spInterfaceHandlerLocal->get_index_shift_tri(); + + size_t numAllDoFs = m_numGridNodes; + + if ( shift_global_index_tri ) + numAllDoFs = m_numGridNodes + m_numNewDoFs; + +// loop all entries for mapping + for (size_t dof1 = 0; dof1 < locJ_tri.num_all_row_dof(0); ++dof1) + { + // get the real index: this data was stored in 'InterfaceHandlerDiffusion::m_vRealCornerID_' + // case1: node lies on interface => real_dof = entry within map 'InterfaceHandlerDiffusion::m_MapNearVertices' + // case2: node lies on an original mesh node: => real_dof = usual, local index of vertex + const size_t dof1_real = m_spInterfaceHandlerLocal->real_index_tri(dof1); + + // case1: dof1_real is index of m_vertex (i.e. interface-DoF): compute global index: + if ( m_spInterfaceHandlerLocal->lies_onInterface_tri(dof1) ) + indexRow = DoFIndex(numAllDoFs + dof1_real, 0); + else // case2: use dof1_real as local index + indexRow = DoFIndex(rowInd.index(0, dof1_real), rowInd.comp(0, dof1_real)); + + for (size_t dof2 = 0; dof2 < locJ_tri.num_all_col_dof(0); ++dof2) + { + const size_t dof2_real = m_spInterfaceHandlerLocal->real_index_tri(dof2); + + // if dof1_real is index of m_vertex (i.e. interface-DoF): compute global index: + if ( m_spInterfaceHandlerLocal->lies_onInterface_tri(dof2) ) + indexCol = DoFIndex(numAllDoFs + dof2_real, 0); + else + indexCol = DoFIndex(colInd.index(0, dof2_real), colInd.comp(0, dof2_real)); + + // finally add loc to glob: + DoFRef(mat, indexRow, indexCol) += locJ_tri.value(0, dof1, 0, dof2); + } + } + +//////////////////////////////////////////////////////////////////////////////////// +/// SECOND: add the boundary contribution stored in 'locJ_quad' to global matrix: +//////////////////////////////////////////////////////////////////////////////////// + + LocalMatrix locJ_quad = get_local_jacobian(ROID_QUADRILATERAL); + const bool shift_global_index_quad = m_spInterfaceHandlerLocal->get_index_shift_quad(); + +// reset numAllDoFs! + numAllDoFs = m_numGridNodes; + + if ( shift_global_index_quad ) + numAllDoFs = m_numGridNodes + m_numNewDoFs; + +// loop all entries for mapping + for (size_t dof1 = 0; dof1 < locJ_quad.num_all_row_dof(0); ++dof1) + { + size_t dof1_real = m_spInterfaceHandlerLocal->real_index_quad(dof1); + + // if dof1_real is index of m_vertex: compute global index: + if ( m_spInterfaceHandlerLocal->lies_onInterface_quad(dof1) + && !m_spInterfaceHandlerLocal->check_vertex_modus(ON_INTERFACE, dof1, -1) ) + indexRow = DoFIndex(numAllDoFs + dof1_real, 0); + else + indexRow = DoFIndex(rowInd.index(0, dof1_real), rowInd.comp(0, dof1_real)); + + for (size_t dof2 = 0; dof2 < locJ_quad.num_all_col_dof(0); ++dof2) + { + size_t dof2_real = m_spInterfaceHandlerLocal->real_index_quad(dof2); + + // if dof1_real is index of m_vertex: compute global index: + if ( m_spInterfaceHandlerLocal->lies_onInterface_quad(dof2) + && !m_spInterfaceHandlerLocal->check_vertex_modus(ON_INTERFACE, dof2, -1) ) + indexCol = DoFIndex(numAllDoFs + dof2_real, 0); + else + indexCol = DoFIndex(colInd.index(0, dof2_real), colInd.comp(0, dof2_real)); + + // finally add loc to glob: + DoFRef(mat, indexRow, indexCol) += locJ_quad.value(0, dof1, 0, dof2); + + } + } + +} + + +template +void DiffusionInterfaceMapper:: +add_local_vec_to_global_Nitsche(vector_type& vec, const LocalVector& lvec, ConstSmartPtr dd) +{ + const LocalIndices& ind = lvec.get_indices(); + DoFIndex index; + + /////////////////////////////////////////////////////////////// + /// FIRST: add loc to glob for 'locD_tri': + /////////////////////////////////////////////////////////////// + LocalVector locD_tri = get_local_defect(ROID_TRIANGLE); + const bool shift_global_index_tri = m_spInterfaceHandlerLocal->get_index_shift_tri(); + + if ( lvec.num_all_dof(0) != locD_tri.num_all_dof(0) ) + UG_THROW("in 'add_local_vec_to_global_Nitsche': non-consistent sizees!\n"); + + size_t numAllDoFs = m_numGridNodes; + + if ( shift_global_index_tri ) + { numAllDoFs = m_numGridNodes + m_numNewDoFs; + } + + for (size_t dof = 0; dof < locD_tri.num_all_dof(0); ++dof) + { + size_t global_index = ind.index(0, dof); + + // if dof_real is index of m_vertex: compute global index: + if ( m_spInterfaceHandlerLocal->lies_onInterface(dof) ) + global_index = numAllDoFs + m_spInterfaceHandlerLocal->get_index_for_global_index_Nitsche(global_index); + + index = DoFIndex(global_index, ind.comp(0, dof)); + + // finally add loc to glob: + DoFRef(vec, index) += locD_tri.value(0, dof); + + } + + + /////////////////////////////////////////////////////////////// + /// SECOND: add loc to glob for 'lvec': + /// -> same as FIRST, but: + /// (1) lvec instead of locD_tri + /// (2) if ( !m_spInterfaceHandlerLocal->lies_onInterface(dof) ) + /// instead of + /// if ( m_spInterfaceHandlerLocal->lies_onInterface(dof) ) + /////////////////////////////////////////////////////////////// + + for (size_t dof = 0; dof < lvec.num_all_dof(0); ++dof) + { + size_t global_index = ind.index(0, dof); + + // if dof_real is index of m_vertex: compute global index: + if ( !m_spInterfaceHandlerLocal->lies_onInterface(dof) ) + global_index = numAllDoFs + m_spInterfaceHandlerLocal->get_index_for_global_index_Nitsche(global_index); + + index = DoFIndex(global_index, ind.comp(0, dof)); + + // finally add loc to glob: + DoFRef(vec, index) += lvec.value(0, dof); + + } + +} + +template +void DiffusionInterfaceMapper:: +add_local_vec_to_global_interface(vector_type& vec, const LocalVector& lvec, ConstSmartPtr dd) +{ + DoFIndex index_print; + + const LocalIndices& ind = lvec.get_indices(); + DoFIndex index; + +/////////////////////////////////////////////////////////////////////////////////// +/// FIRST: add the boundary contribution stored in 'locD_tri' to global defect: +/////////////////////////////////////////////////////////////////////////////////// + + LocalVector locD_tri = get_local_defect(ROID_TRIANGLE); + const bool shift_global_index_tri = m_spInterfaceHandlerLocal->get_index_shift_tri(); + + size_t numAllDoFs = m_numGridNodes; + + if ( shift_global_index_tri ) + numAllDoFs = m_numGridNodes + m_numNewDoFs; + +// UG_LOG("in DiffusionInterfaceMapper::add_loc_vec(): locD_tri = " << locD_tri << "\n"); + + for (size_t dof = 0; dof < locD_tri.num_all_dof(0); ++dof) + { + const size_t dof_real = m_spInterfaceHandlerLocal->real_index_tri(dof); + + // if dof_real is index of m_vertex: compute global index: + if ( m_spInterfaceHandlerLocal->lies_onInterface_tri(dof) && !m_spInterfaceHandlerLocal->check_vertex_modus(ON_INTERFACE, dof, 1) ) + index = DoFIndex(numAllDoFs + dof_real,0); + else + index = DoFIndex(ind.index(0, dof_real), ind.comp(0, dof_real)); + + // finally add loc to glob: + DoFRef(vec, index) += locD_tri.value(0, dof); + + } + + +///////////////////////////////////////////////////////////////////////////////// +/// SECOND: add the boundary contribution stored in 'locD_quad' to global defect: +///////////////////////////////////////////////////////////////////////////////// + + LocalVector locD_quad = get_local_defect(ROID_QUADRILATERAL); + const bool shift_global_index_quad = m_spInterfaceHandlerLocal->get_index_shift_quad(); + +// reset numAllDoFs! + numAllDoFs = m_numGridNodes; + + if ( shift_global_index_quad ) + numAllDoFs = m_numGridNodes + m_numNewDoFs; + +// UG_LOG("in DiffusionInterfaceMapper::add_loc_vec(): locD_quad = " << locD_quad << "\n"); + + for (size_t dof = 0; dof < locD_quad.num_all_dof(0); ++dof) + { + size_t dof_real = m_spInterfaceHandlerLocal->real_index_quad(dof); + + // if dof_real is index of m_vertex: compute global index: + if ( m_spInterfaceHandlerLocal->lies_onInterface_quad(dof) + && !m_spInterfaceHandlerLocal->check_vertex_modus(ON_INTERFACE, dof, 1) ) + index = DoFIndex(numAllDoFs + dof_real,0); + else + index = DoFIndex(ind.index(0, dof_real), ind.comp(0, dof_real)); + + // finally add loc to glob: + DoFRef(vec, index) += locD_quad.value(0, dof); + } + + +} + +} // namespace ConvectionDiffusionPlugin +} // end ug namespace + +#endif /* DIFFUSION_INTERFACE_MAPPER_IMPL_H_ */