Sanne conv diff potential master

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

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<TDomain> T;
+        typedef ConvectionDiffusionBase<TDomain> TBase;
+        string name = string("ConvectionDiffusionFV1_cutElem").append(suffix);
+        reg.add_class_<T, TBase >(name, grp)
+        .template add_constructor<void (*)(const char*,const char*)>("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<TDomain> 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 <typename TDomain, typename TAlgebra>
+    static void DomainAlgebra(Registry& reg, string grp)
+    {
+        string suffix = GetDomainAlgebraSuffix<TDomain,TAlgebra>();
+        string tag = GetDomainAlgebraTag<TDomain,TAlgebra>();
+
+        typedef ApproximationSpace<TDomain> approximation_space_type;
+        typedef GridFunction<TDomain, TAlgebra> function_type;
+
+
+        //	ImmersedInterfaceDiffusion
+        {
+
+            typedef ImmersedInterfaceDiffusion<TDomain, TAlgebra> T;
+            typedef IImmersedInterface<TDomain, TAlgebra> TBase;
+            string name = string("ImmersedInterfaceDiffusion").append(suffix);
+            reg.add_class_<T, TBase>(name, grp)
+            .template add_constructor<void (*)(SmartPtr<IAssemble<TAlgebra> > ass,
+                                               SmartPtr<ConvectionDiffusionPlugin::ConvectionDiffusionFV1_cutElem<TDomain> > spMaster,
+                                               SmartPtr<DiffusionInterfaceProvider<TDomain::dim> > interfaceProvider,
+                                               SmartPtr<CutElementHandler_TwoSided<TDomain::dim> > 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 <int dim>
 static void Dimension(Registry& reg, string grp)
 {
 	string dimSuffix = GetDimensionSuffix<dim>();
 	string dimTag = GetDimensionTag<dim>();
-
-	//	singular sources and sinks 
+    
+	//	singular sources and sinks
 	{
 		typedef CDSingularSourcesAndSinks<dim> 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<dim> T;
+        string name = string("CutElementHandler_TwoSided").append(suffix);
+        reg.add_class_<T>(name, grp)
+        .template add_constructor<void (*)(SmartPtr<MultiGrid> mg, const char*, SmartPtr<DiffusionInterfaceProvider<dim> >)>("multigrid, fct names")
+        .set_construct_as_smart_pointer(true);
+        reg.add_class_to_group(name, "CutElementHandler_TwoSided", tag);
+    }
+
+    // DiffusionInterfaceProvider
+    {
+        typedef DiffusionInterfaceProvider<dim> T;
+        string name = string("DiffusionInterfaceProvider").append(suffix);
+        reg.add_class_<T>(name, grp)
+        .template add_constructor<void (*)( )>("")
+        .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<Functionality>(*reg,grp);
 		RegisterDomainDependent<Functionality>(*reg,grp);
+        RegisterDomainAlgebraDependent<Functionality>(*reg,grp);
 		RegisterDomain2d3dDependent<Functionality2d3d>(*reg,grp);
 	}
 	UG_REGISTRY_CATCH_THROW(grp);

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

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.
+
+
+