GlobalBasis.md 13.3 KB
 Müller, Felix committed Sep 29, 2020 1 2 # Global Basis {: #group-globalbasis } ## Summary  Praetorius, Simon committed Nov 07, 2020 3 In the context of the finite-element method we use a finite-element space $V$ with the set of basis functions $\{\phi_i\}$. Within AMDiS this concept is realized by the class GlobalBasis, contained in the header file amdis/functions/GlobalBasis.hpp. This extends the interface of the underlying Dune::Functions::DefaultGlobalBasis with an automatic update mechanism used in several places within AMDiS. We strongly advice to always use an AMDiS::GlobalBasis or a user-defined derived class instead of the DUNE data structure for the update mechanism to work properly.  Müller, Felix committed Sep 29, 2020 4 5 6 7 8 For more information on the class interface visit the API documentation. ## PreBasis classes and basis trees Many finite-element spaces in applications can be constructed as a product of simple spaces $V = V_1 \times V_2 \times \dots \times V_k$. For example the Taylor-Hood-Element can be constructed as $V_{TH} = V_{v} \times V_{p}$ with the space of velocity functions $V_v$ and pressure functions $V_p$. The velocity space can again be decomposed into the vector components $V_v = V_{v_1} \times \dots \times V_{v_n}$. If we use second-degree lagrange basis functions for the velocity space and first-order lagrange basis functions for the pressure we get a decomposition $V_{TH} = V_{v_1} \times V_{v_n} \times V_p = L_2^n \times L_1$.  Praetorius, Simon committed Nov 07, 2020 9 The underlying numerics environment of AMDiS, DUNE, attempts to model the situation above. Hence a GlobalBasis in AMDiS can be defined by a tree structure using the Dune-PreBasis class, which itself is composed of nodes in a tree. The leaf nodes in this tree contain the implementation details of the simple basis, while inner nodes indicate composition of either identical children (called power nodes) or arbitrary children (called composite nodes).  Müller, Felix committed Sep 29, 2020 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39  ### Making a PreBasis When we want to build a PreBasis we must first identify the structure of the basis we want to construct with a tree. We can then build that structure in AMDiS by nesting the functions Dune::Functions::BasisFactory::composite(Nodes...[, MergingStrategy]) for composite nodes, Dune::Functions::BasisFactory::power(Node[, MergingStrategy]) for power nodes and implementations (e.g. Dune::Functions::BasisFactory::lagrange())). The second optional argument MergingStrategy provides merging strategies to the inner nodes, specifying how the indices of the simple leaf basis should be merged to obtain the indexing of the global basis. Currently only flat indexing is supported by AMDiS. The following code snippet shows how a PreBasisFactory for a Taylor-Hood-Element is constructed, that can later be used to build a global basis. c++ using namespace Dune::Functions::BasisFactory; const int dow = 2; // world dimension const int k = 1; // order parameter auto taylorHoodPreBasisFactory = composite( power( lagrange() ), lagrange() ); auto taylorHoodPreBasisFactoryWithMergingStrategy = composite( power( lagrange(), flatInterleaved() ), lagrange(), flatLexicographic() );  ## Making a global basis  Praetorius, Simon committed Nov 07, 2020 40 Using a PreBasisFactory we can easily make a global basis by defining the set of grid elements on which the basis functions of the FE-space should live. This can be done by providing a GridView and using the GlobalBasis constructors. An optional name can be provided that can be used to pass initfile parameters to the parallel communication class.  Müller, Felix committed Sep 29, 2020 41 42 43 44  c++ // Name, Grid and Dune::DefaultGlobalBasis arguments template  Praetorius, Simon committed Nov 07, 2020 45  GlobalBasis(std::string const& name, Grid const& grid, Args&&... args);  Müller, Felix committed Sep 29, 2020 46 47 48  // As above with name defaulting to "" template  Praetorius, Simon committed Nov 07, 2020 49  GlobalBasis(Grid const& grid, Args&&... args);  Müller, Felix committed Sep 29, 2020 50 51 52  // Name, GridView and PreBasisFactory template  Praetorius, Simon committed Nov 07, 2020 53  GlobalBasis(std::string const& name, GridView const& gridView, PBF const& preBasisFactory);  Müller, Felix committed Sep 29, 2020 54 55 56  // As above with name defaulting to "" template  Praetorius, Simon committed Nov 07, 2020 57  GlobalBasis(GridView const& gridView, PBF const& preBasisFactory);  Müller, Felix committed Sep 29, 2020 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85  If we use a ProblemStat object most of the work will be done automatically for us. The PreBasis is specified via a Traits template parameter with the most frequently used cases already included in amdis/ProblemStatTraits.hpp. Those include c++ // Composition of any number of lagrange bases with any degree template struct LagrangeBasis; // As above but a structured grid is chosen template struct YaspGridBasis; // The taylor-hood basis as discussed above template struct TaylorHoodBasis;  If one of the above traits class is provided ProblemStat will create a new global basis on a call to ProblemStat::initialize. Afterwards a pointer to the basis can be obtained using the function ProblemStat::globalBasis. c++ using Grid = Dune::YaspGrid<2>; // 2-dimensional structured grid ProblemStat> prob("myProblem"); prob.initialize(INIT_ALL); auto& basis = *prob.globalBasis();  ## Using the global basis  Praetorius, Simon committed Nov 07, 2020 86 The GlobalBasis provides access to element indices and basis functions in the same way as a Dune-basis does, it is even derived from the Dune::DefaultGlobalBasis.  Müller, Felix committed Sep 29, 2020 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124  ### Getting the total number of DOFs If we are simple interested in getting the total number of DOFs the basis contains we can simply call basis.dimension(). This can be useful for preparing matrices and vectors for insertion. c++ std::size_t maxSize = basis.dimension(); std::vector dofs; dofs.resize(maxSize);  ### Access to local basis functions using a LocalView Within that interface we are restricted to elementwise access using a LocalView. This provides us with a way to work with the local basis functions on one grid element. We show a typical use in the following snippet. Note how we first need to bind the LocalView before we can use it. c++ auto localView = basis.localView(); for (const auto& e : elements(basis.gridView())) // loop over all grid elements the basis is defined on { localView.bind(e); // a LocalView must be bound to an element before being used // do something localView.unbind(); }  A bound LocalView has the method LocalView::index(size_type) mapping a local index to a global index. In other words it maps a local basis function defined on an element to its corresponding global basis function. We can use that to build a global stiffness matrix from local contributions on a single element and then insert those into a single matrix in global indices. Another method is LocalView::tree() that returns the root node of the local basis tree. The main method all nodes share is Node::localIndex(size_type) which maps a leaf node index to the local index within the local basis tree. #### The for_each_node and for_each_leaf_node helper functions Quite often we want to perform operations on certain nodes of the tree other than the root node. This can be useful if we want to work with the actual implementations wich are usually leaf nodes. For this we can use the helper functions for_each_node and for_each_leaf_node defined in amdis/typetree/Traversal.hpp. Those functions traverse the tree and call the given function on every (leaf) node with the node and a type of tree index we shall explain later as arguments. we show the usage with the following example using the Taylor-Hood-Basis defined above. Here we assume to have a LocalView localView that is bound to an element. c++ auto root = localView.tree(); for_each_leaf_node(root, [&](auto const& node, auto const& tp) { // do something on node });  #### Working on specific nodes using a TreePath {: #globalbasis-using-treepath } There are cases when we want to address a certain tree node. To come back to our Taylor-Hood example we might want to add an operator that only acts on the velocity vector or the pressure, or even just a single component of the velocity vector. For this end a type of node index exists, called TreePath. This defines a list of indices specifying which path to take from the root at each node on the way. For technical reasons the index type is either an int (for power nodes) or std::integral_constant for composite nodes. For convenience we can use Dune::Indices::_0 instead of writing std::integral_constant. Note how indices always start at 0 and follow the order we specify when creating the PreBasis.  Praetorius, Simon committed Oct 14, 2020 125 Using once again the Taylor-Hood basis we can use the makeTreePath(Args...) function to convert indices into a TreePath for certain nodes.  Müller, Felix committed Sep 29, 2020 126 127  c++  Praetorius, Simon committed Oct 14, 2020 128 129 130 131 auto tp_v = makeTreePath(Dune::Indices::_0); // velocity vector auto tp_p = makeTreePath(Dune::Indices::_1); // pressure auto tp_v2 = makeTreePath(Dune::Indices::_0, 1); // second component of the velocity vector auto tp_root = makeTreePath(); // root  Müller, Felix committed Sep 29, 2020 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200  Recall that localView.tree() returns the root of the basis tree. Using a TreePath we can access a specific node using the free funtion Dune::TypeTree::child(Node, TreePath). c++ auto root = localView.tree(); auto node_v = Dune::TypeTree::child(root, tp_v); // power node representing the velocity vector auto node_p = Dune::TypeTree::child(root, tp_p); // lagrange node representing the pressure  ### Using the lagrange (pre-)basis Many applications require only lagrange elements or compositions, for example a Navier-Stokes problem may use the Taylor-Hood basis we introduced above that consists of dim second order lagrange elements and one first order lagrange element. For that reason we will now take a closer look at the implementation of those lagrange elements in AMDiS. AMDiS borrows the implementation from the underlying DUNE-functions module. This defines lagrange elements of a given order k. Recall that we can add lagrange nodes to a (pre-)basis using Dune::Functions::BasisFactory::lagrange(). In the previous section we have seen how we can access the leaf nodes of a basis. With that we can get the implementation of the local finite element with the function Dune::Functions::LagrangeNode::finiteElement(). We shall show the usage of the local finite element class handed out by the function above. Its interface defines the functions size(), localCoefficients(), localInterpolation() and localBasis(). We shall explain those with an example. Assume we have localView bound to an element and have vectors dofs and dofs2 that store the coefficients of some grid function with global indexing (e.g. a solution vector). Domain is the type used for local coordinates of the reference element. c++ auto root = localView.tree(); for_each_leaf_node(root, [&](auto const& node, auto const& tp) { // Extract some types from the node using Node = Underlying_t; using LocalFunction = typename Node::FiniteElement::Traits::LocalInterpolationType::FunctionType; using Domain = typename LocalFunction::Traits::DomainType; using Range = typename LocalFunction::Traits::RangeType; auto const& fe = node.finiteElement(); auto feSize = fe.size(); // get information about the position of a basis function in the element auto const& localCoefficients = fe.localCoefficients(); for (std::size_t i = 0; i < feSize; ++i) { auto const& localKey = localCoefficients.localKey(i); // do something } // interpolate a local function onto an element auto f_ = [](Domain const& x) -> Range { return x[0] + x[1]; }; auto f = functionFromCallable(f_); std::vector localCoeff; fe.localInterpolation().interpolate(f, localCoeff); // interpolate f onto the local basis for (std::size_t i = 0; i < feSize; ++i) { std::size_t globalIndex = localView.index(node.localIndex(i)); // get the global index dofs[globalIndex] = localCoeff[i]; // set global coefficient to the local coefficient } // evaluate at a local coordinate std::vector localContrib(feSize); fe.localBasis().evaluateFunction({0.2, 0.2}, localContrib); // get contributions of local basis functions at a given local coordinate Range sum = 0; for (std::size_t i = 0; i < feSize; ++i) { std::size_t globalIndex = localView.index(node.localIndex(i)); // get the global index sum += dofs2[globalIndex] * localContrib[i]; // add contribution from i-th local basis function } });  ### Keeping indices and data updated when the grid changes Within an adaptive simulation we may want to add or remove grid elements by refinement or coarsening. When this happens the number of total elements or their relative position may change. Therefore the indexing scheme used by the global basis must be updated when changes to the underlying grid happen. For this purpose the update method exists, which takes a GridView of the Grid after it has been changed as its only argument. c++ void update(GridView const& gv)  Usually we do not need to call this function - out of the box AMDiS with a ProblemStat doing most of the work will automatically call it for us. If user code is working with the underlying Dune-Grid directly there is no way for AMDiS to detect if any changes happen to it, therefore update must be called manually in such a scenario. ### The globalRefineCallback function  Praetorius, Simon committed Nov 07, 2020 201 Certain grid managers support the use of a callback function when doing global refinement. Using a GlobalBasis in this context is currently not supported.