Sanne naver stokes potential master

CMakeLists.txt

@@ -52,7 +52,8 @@ set(SOURCES	navier_stokes_base.cpp
 			incompressible/incompressible_navier_stokes_base.cpp	
 
 			incompressible/fv1/navier_stokes_fv1.cpp
-			incompressible/fe/navier_stokes_fe.cpp
+            incompressible/fv1/navier_stokes_fv1_cutElem.cpp
+            incompressible/fe/navier_stokes_fe.cpp
 			incompressible/fvcr/navier_stokes_fvcr.cpp
 			incompressible/fv/navier_stokes_fv.cpp
 
@@ -74,7 +75,8 @@ set(SOURCES	navier_stokes_base.cpp
 			incompressible/fvcr/register_fvcr.cpp
 			incompressible/fe/register_fe.cpp	
 
-			incompressible/incompressible_navier_stokes_plugin.cpp)
+			incompressible/incompressible_navier_stokes_plugin.cpp
+            incompressible/particle_laden_flow_plugin.cpp)
 
 
 ################################################################################

incompressible/CMakeLists.txt

@@ -51,6 +51,7 @@ set(SOURCES	../navier_stokes_base.cpp
 			incompressible_navier_stokes_base.cpp	
 
 			fv1/navier_stokes_fv1.cpp
+            fv1/navier_stokes_fv1_cutElem.cpp
 			fe/navier_stokes_fe.cpp
 			fvcr/navier_stokes_fvcr.cpp
 			fv/navier_stokes_fv.cpp

incompressible/fv1/bnd/inflow_fv1_cutElem.h

@@ -0,0 +1,101 @@
+/*
+ * Copyright (c) 2013:  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__NAVIER_STOKES__INCOMPRESSIBLE__FV1__BND__INFLOW_FV1_CUT_ELEM__
+#define __H__UG__NAVIER_STOKES__INCOMPRESSIBLE__FV1__BND__INFLOW_FV1_CUT_ELEM__
+
+#include "../../bnd/inflow_base.h"
+#include "../navier_stokes_fv1_cutElem.h"
+
+#include "lib_disc/spatial_disc/elem_disc/neumann_boundary/neumann_boundary_base.h"
+#include "lib_disc/spatial_disc/constraints/dirichlet_boundary/lagrange_dirichlet_boundary.h"
+#include "lib_disc/spatial_disc/elem_disc/neumann_boundary/fv1/neumann_boundary_fv1.h"
+
+namespace ug{
+namespace NavierStokes{
+
+template <	typename TDomain, typename TAlgebra>
+class NavierStokesInflowFV1_cutElem
+	: public NavierStokesInflowBase<TDomain, TAlgebra>
+{
+	private:
+		static const int dim = TDomain::dim;
+
+	public:
+
+	///	returns the number of element discs
+		virtual size_t num_elem_disc() const {return 1;}
+
+	///	returns the element disc
+		virtual SmartPtr<IElemDisc<TDomain> > elem_disc(size_t i) {return m_spNeumannDisc;}
+
+	///	returns the number of constraints
+		virtual size_t num_constraint() const {return 1;}
+
+	///	returns an element disc
+		virtual SmartPtr<IDomainConstraint<TDomain, TAlgebra> > constraint(size_t i) {return m_spDirichletConstraint;}
+
+	public:
+	///	Constructor
+    NavierStokesInflowFV1_cutElem(SmartPtr< NavierStokesFV1_cutElem<TDomain> > spMaster)
+    : m_spNeumannDisc(new NeumannBoundaryFV1<TDomain>(spMaster->symb_fcts()[dim].c_str())),
+    m_spDirichletConstraint(new DirichletBoundary<TDomain,TAlgebra>),
+    m_spMaster(spMaster) // copy constructor not defined for _cutElem?
+    {
+        const std::vector<std::string>& vFctName = m_spMaster->symb_fcts();
+
+        if(vFctName.size() != TDomain::dim + 1)
+            UG_THROW("NavierStokesInflow::set_functions: This Boundary "
+                     "Condition works on exactly dim+1 (velocity+pressure) "
+                     "components, but "<<vFctName.size()<<"components given.");
+    }
+
+	///	sets the velocity to a given value
+		void add(SmartPtr<CplUserData<MathVector<dim>, dim> > user, const char* subsetsBND);
+
+	protected:
+	///	neumann disc for pressure equation
+		SmartPtr<NeumannBoundaryBase<TDomain> > m_spNeumannDisc;
+
+	///	dirichlet disc for velocity components
+		SmartPtr<DirichletBoundary<TDomain,TAlgebra> > m_spDirichletConstraint;
+
+	/// The master discretization:
+		SmartPtr< NavierStokesFV1_cutElem<TDomain> > m_spMaster;
+};
+
+} // namespace NavierStokes
+} // end namespace ug
+
+#include "inflow_fv1_cutElem_impl.h"
+
+#endif /* __H__UG__NAVIER_STOKES__INCOMPRESSIBLE__FV1__BND__INFLOW_FV1_CUT_ELEM__ */

incompressible/fv1/bnd/inflow_fv1_cutElem_impl.h

@@ -0,0 +1,87 @@
+/*
+ * Copyright (c) 2013:  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__NAVIER_STOKES__BND__INFLOW_FV1_CUT_ELEM_IMPL__
+#define __H__UG__NAVIER_STOKES__BND__INFLOW_FV1_CUT_ELEM_IMPL__
+
+#include "inflow_fv1_cutElem.h"
+#include "lib_disc/spatial_disc/elem_disc/neumann_boundary/fv1/neumann_boundary_fv1.h"
+
+namespace ug{
+namespace NavierStokes{
+/*
+template <typename TDomain, typename TAlgebra>
+NavierStokesInflowFV1_cutElem<TDomain,TAlgebra>::
+NavierStokesInflowFV1_cutElem(SmartPtr< NavierStokesFV1_cutElem<TDomain> > spMaster)
+	: m_spNeumannDisc(new NeumannBoundaryFV1<TDomain>(spMaster->symb_fcts()[dim].c_str())),
+	  m_spDirichletConstraint(new DirichletBoundary<TDomain,TAlgebra>),
+	  m_spMaster(spMaster)
+{
+	const std::vector<std::string>& vFctName = m_spMaster->symb_fcts();
+
+	if(vFctName.size() != TDomain::dim + 1)
+		UG_THROW("NavierStokesInflow::set_functions: This Boundary "
+				"Condition works on exactly dim+1 (velocity+pressure) "
+				"components, but "<<vFctName.size()<<"components given.");
+}
+*/
+template <typename TDomain, typename TAlgebra>
+void NavierStokesInflowFV1_cutElem<TDomain,TAlgebra>::
+add(SmartPtr<CplUserData<MathVector<dim>, dim> > user, const char* subsetsBND)
+{
+	const std::vector<std::string>& vFctName = m_spMaster->symb_fcts();
+	if(vFctName.empty())
+		UG_THROW("NavierStokesInflow::add: Symbolic names for"
+				" velocity and pressure not set. Please set them first.");
+
+	std::string innerSubsets;
+	for(size_t s = 0; s < m_spMaster->symb_subsets().size(); ++s){
+		if(s > 0) innerSubsets.append(",");
+		innerSubsets.append(m_spMaster->symb_subsets()[s]);
+	}
+
+
+	m_spNeumannDisc->add(user, subsetsBND, innerSubsets.c_str());
+
+	std::string velNames;
+	for(int i=0;i<TDomain::dim; ++i)
+	{
+		if(i>0) velNames.append(",");
+		velNames.append(vFctName[i]);
+	}
+	m_spDirichletConstraint->add(user.template cast_dynamic<UserData<MathVector<dim>, dim> >(), velNames.c_str(), subsetsBND);
+}
+
+} // namespace NavierStokes
+} // end namespace ug
+
+#endif /* __H__UG__NAVIER_STOKES__BND__INFLOW_FV1_CUT_ELEM_IMPL__ */

incompressible/fv1/moving_particle/Info File 1

@@ -0,0 +1,40 @@
+/* Info File for 'MovingParticle' class: */
+
+The 'MovingParticle' class enables the simulation of moving particles in a fluid. Detailed descriptions of the mathematical theory and numerical simulations can be found in [1] and [2].
+The interaction between the fluid and the particles takes place at the interface of the particles, which acts therefore as an immersed interface for the fluid. Consequently, the 'MovingParticle' class implements a derivation of the 'ImmersedInterface' class. It therefore contains all the components of an 'ImmersedInterface' class implemented according to the necessities for moving particles. Those are explained in more detail in the according class types.
+
+
+
+The special requirements for a moving particle immersed interface implementation are the following:
+
+(1) Each particle is described as a rigid motion, i.e. by its two components: the translational and rotational velocity. As a consequence, in the discrete scheme each particle is described by the DoFs for its translational and rotational velocities. Those can be stored in two geometrical nodes of the mesh. Their location within each particle gets defined during an initialisation phase. Most of the nodes within the particles are therefore NO DoFs for the particle, beside the designated nodes (named 'transInd' and 'rotInd' for the associated global indices).
+
+(2) Furthermore, a particle usually covers more than one (local) finite element. Since the solution for the particle velocities is coupled with and need to be known for all surrounding fluid DoFs, the DoFs need to be accessible on a GLOBAL scale as well. The assigned global indices ('transInd' and 'rotInd') for the particle DoFs need to be accessible globally. Furthermore, the initialization of the interfaces of the particles inherits the collection of all cut elements.
+
+
+The two central methods called during the initialisation (called via the 'init()' of the associated 'CutElementHandler') are:
+
+---> update_interface_data(): computes all cut elements
+---> update_global_indices(): defines the nodes and indices within a particle, which act as particle DoFs
+
+
+Further remarks:
+
+- Most of the nodes within the particles are NO DoFs for the particle.
+
+- During the movement of a particle in time, nodes that were part of the particle can get part of the fluid in the next time step. Since they were no designated DoFs within the particle they similarly did not contain the solution. A 'filling' of these so called 'freed' nodes need to be performed during the (geometrical) update of the particle interfaces.
+
+- During the local assembling of the defect, the solution for the translational or rotational components of the particle velocity will be overwritten with the physical particle velocity (a combination of boths). Before and after these computations the particle solution needs to be transfered back and forth between local and global data structures for all later assembling.
+
+
+
+
+[1] Hoellbacher S., Wittum G.:
+    "Rotational test spaces for a fully-implicit FVM and FEM for the DNS of fluid-particle interaction"
+    Journal of Computational Physics, vol. 393, (2019), pp. 186–213
+    open access link: https://authors.elsevier.com/sd/article/S0021999119303298
+    DOI: https://doi.org/10.1016/j.jcp.2019.05.004
+
+[2] Hoellbacher S., Wittum G.:
+    "Gradient-consistent enrichment of finite element spaces for the DNS of fluid-particle interaction."
+    submitted to Jounal of Computational Physics

incompressible/fv1/moving_particle/Info File 2

@@ -0,0 +1,42 @@
+/* Info File for 'ParticleMapper' class: */
+
+The 'ParticleMapper' class handles the mapping from the locally computed data (stiffnes matrix, defect) to the global algebra.
+
+The 2 main components are:
+
+(1) Mapping of the local indices, which are associated to the particle DoFs to their global indices.
+(2) Performing the cross product multiplication for local indices associated to the rotational velocity DoFs.
+
+
+Important methods:
+
+---> the hierarchy
+
+            'add_local_mat_to_global()' ---> 'add_local_mat_to_global_interface()' ---> 'add_local_mat_to_global_FT()'
+
+     forwards the adapted mapping for the cut elements
+
+---> modify_LocalData(): resizes the local data (LocalMatrix or LocalVecor) BEFORE starting the assembling on the element. It is called during the elements discretisation loop (see elem_disc_assemble_util()).
+
+---> modify_GlobalSol(): writes the particle velocities stored in its designated global indices to own data in order to protect it from being overwritting in the following local assembling process and provide global access.
+
+
+
+
+The flat-top 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: Each coupling is mapped to its own, designated global index, i.e. DoF. The resulting space will be the usual finite element space w.r.t. the interface adaped grid. It therefore contains additional nodes along the interface and (virtually) consists of cut elements and original elements.
+
+Case 2: The local coupling 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 same element. The resulting finite element space will potentially be smaller than the one of the cut element, since two and more (in 3d) 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 hyperplane 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. Therefore, the 'ParticleMapper' is a flat-top mapper.
+
+In contrast, the 'DiffusionInterfaceMapper' is of type case 1:
+In particular, it is a tow-sided mapper and the resulting space is an interface adapted space.
+
+
+

incompressible/fv1/moving_particle/Info File 3

@@ -0,0 +1,6 @@
+/* Info File for 'ParticleBndCond' class: */
+
+Assemlbes the stresses on the particle interface. For that, the common forces will be assembled as for usual boundary faces, BUT on the according boundary faces which were computed for the cut element during the local update of the element.
+
+The boundary faces of the cut elements got stored into the 'm_vBF' member of the 'InterfaceHandlerLocalParticle' class.
+