# Session 2 ## Wednesday Nov 28 - Scalar linear second order PDEs - **Discrete functions on unstructured grids** --- # Motivation ### Extract maxima of solution <img src="images/ball500.png" width="25%" /> <img src="images/pfc_sphere.png" width="25%" /> ### Calculate integrals <img src="images/energy.png" width="35%" /> --- # Working with discrete functions Finite-Element function `\(u_h\)`, a linear combination of basis-functions `\(\phi_i\)`: \\[ u\_h(x) = \sum\_i u\_i \phi\_i(x) = \sum\_j u\_{i\_T(j)} \hat{\phi}_j(\lambda\_T(x)), \\] with `\(u_i:=\)` Degrees of Freedom, stored in `DOFVector<value_type> U`: - Evaluate at DOF index: `u_i = U[i]` - Iterate over all DOFs: ``` DOFIterator<value_type> it(&U [, USED_DOFS]); for (it.reset(); !it.end(); ++it) *it = value; // idx = it.getDOFIndex() ``` - Evaluate DOFVector at global coordinate `\(x\)`: ``` WorldVector<double> x; x[0] = 0.2; x[1] = 0.3; value_type value = U(x); ``` --- ## Discrete function on element T Evaluate DOFVector in local coordinates `\(\lambda\)` of element `T=elInfo->getElement()` \\[ u\_h(x\_T(\lambda)) = \sum\_j \hat{u}\_j \hat{\phi}_j(\lambda), \\] with `\(\hat{u}_j = u_{i_T(j)}\)` 1. Get local DOFs `u_j` on the element 2. Evaluate the linear combination of basis functions ``` const FiniteElemeSpace* feSpace = U.getFeSpace(); const BasisFunction* basFcts = feSpace->getBasisFcts(); ElementVector uh(basFcts->getNumber()); U.getLocalVector(elInfo->getElement(), uh); // local DOFs u_j basFcts->evalUh(lambda, uh); // linear combination ``` --- # Interpolate a local function to a DOFVector A local function is something that can be evaluated at local coordinates in a *bound* element. By traversing the mesh the local function can be used to locally interpolate to the DOFs of an element. ``` U << EXPRESSION; ``` where `EXPRESSION` can contain e.g. - The coordinates: `X()`, `X(i)` - Scalar values: `1.0`, `constant(1.0)` - Other DOFVectors: `valueOf(V)`, `gradientOf(V)` - Matrix and vector expressions: `two_norm(...)`, `vec * vec` ([Expressions manual](https://gitlab.mn.tu-dresden.de/iwr/amdis/wikis/expressions)) --- ### Example: Let `\(C\)` and `\(U\)` be of type `DOFVector<double>`: \\[ U(x) := \frac{1}{2}\left(C(x)^2\,(1-C(x))^2 + \frac{1}{\epsilon}\|\nabla C(x)\|^2\right) \\] Assign expression on rhs to DOFVector `\(U\)`: ``` U << 0.5*( pow<2>(valueOf(C) * (1.0 - valueOf(C)) + (1.0/eps) * unary_dot(gradientOf(C)) ); ``` --- # Working with the discrete solution - Reduce the `DOFVector`, i.e. calc integrals, norms: ``` integrate( EXPRESSION [, Mesh*] ) ``` ``` integrate( valueOf(U) ); integrate( two_norm(gradientOf(U)) ); integrate(pow<2>(valueOf(U)) + unary_dot(gradientOf(U))); ``` - Other reduction operations (over DOFs): - `max( EXPRESSION )` `\(\Rightarrow\max\{ expr(x_i)\,:\, x_i\in\Omega \}\)` - `min( EXPRESSION )` `\(\Rightarrow\min\{ expr(x_i)\,:\, x_i\in\Omega \}\)` - `abs_max( EXPRESSION )` `\(\Rightarrow\max\{ |expr(x_i)|\,:\, x_i\in\Omega \}\)` - `abs_min( EXPRESSION )` `\(\Rightarrow\min\{ |expr(x_i)|\,:\, x_i\in\Omega \}\)` --- # Working with the mesh .center[ <img src="images/ellipt_macro.png" width="30%" /><img src="images/torus_macro.jpg" width="30%" /> ] ### Change position of mesh vertices: <img src="images/Fsi.png" width="70%" /> --- # Working with the mesh A Mesh holds all information about a triangulation. Get the mesh from `ProblemStat`: ``` Mesh& mesh = *prob.getMesh(); ``` - Number of vertices: `mesh.getNumberOfVertices()` - Number of triangles (elements): `mesh.getNumberOfElements()` - World coordinates are stored in `WorldVector<double>`. - Get the coordinates of all DOFs: ``` DOFVector<WorldVector<double> > Coords(prob.getFeSpace(), "coords"); mesh.getDofCoords(Coords); ``` - Change coordinates of the mesh: ``` ParametricSimple parametric(Coords); mesh.setParametric(&parametric); ``` --- # Working with the mesh: - (Advanced) Traverse the mesh element-wise: ``` Flag traverseFlag = Mesh::CALL_LEAF_EL | (FILL_FLAGS); TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(&mesh, -1, traverseFlag); while (elInfo) { // do something with elInfo elInfo = stack.traverseNext(elInfo); } ``` - (Advanced) Global Refinement of the mesh: Use `RefinementManager` ``` RefinementManager* refManager = prob.getRefinementManager(); Flag f = refManager->globalRefine(&mesh, 5); if (f == MESH_REFINED) { MSG("Mesh globally refined by 5 levels!\n"); } ``` --- # Example: Print solution with coordinates Writing out the solution to the screen: ``` DOFVector<WorldVector<double>> C(U.getFeSpace(), "c"); mesh.getDofCoords(C); DOFIterator<double> it(U); DOFIterator<WorldVector<double>> c_it(C); for (it.reset(),c_it.reset(); !it.end(); ++it,++c_it) std::cout << "U(" << *c_it << ") = " << *it << "\n"; ``` This will print out: <pre><code style="background: #ddd;">U(0 0) = 0 U(1 0) = 0 U(1 1) = 0 U(0 1) = 0 U(0.5 0.5) = 0.0833333 </code></pre> --- class: center, middle # Input/Output --- # Input/Output Write DOFVectors to file, for visualization, serialization, ... - File writer with automatic file-type detection: ``` io::writeFile(DOFVECTOR, FILENAME); ``` Detection by filename extension: ".vtu", ".arh", ".dat", ".2d", ".gnu",... - Use **ParaView** for visualization - Tool **MeshConv** can convert beween mesh formats - File reader with automatic file-type detection: ``` io::readFile(FILENAME, DOFVECTOR); ``` - Use ARH format to exchange data (can be converted to VTU by MeshConv.) --- class: center, middle # Exercise2 ## Poisson equation --- # Exercise2: Poisson equation We want to solve the Poisson equation \\[ -\Delta u = f(x)\quad\text{ in }\Omega,\quad u|\_{\partial\Omega} = g \\] in a rectangular domain `\(\Omega\)`, with \\[ f(x) = 40(1 - 10\|x\|^2)e^{-10.0\|x\|^2}\\\\ g(x) = e^{-10.0\|x\|^2} \\\\ \\] and with exact solution `\(u^\ast=g\)`. 1. Assemble and solve the equation. 2. Interpolate error `\(ERR := |u_h - g|\)` to DOFVector and write it to a file. 3. Calculate the error-norms `\(\|u_h - g\|_{L_2(\Omega)}\)` and `\(\|u_h - g\|_{H_1(\Omega)}\)` 4. Evaluate the error at the grid-point `\(x=(0.5, 0.5)\)`. 5. Refine the mesh globally and compare error-norms to old error-norms. --- # Advanced Exercise2: Mesh adaption The datastructure `AdaptInfo` provides parameters for tolerance and nr. of iterations for an adaption process: ``` AdaptInfo adaptInfo("adapt"); adaptInfo.getMaxSpaceIteration(); // => adapt->max iteration adaptInfo.getSpaceTolerance(0); // => adapt[0]->tolerance ``` 1. Use AdaptInfo to write an adaption loop that refines the mesh globally to reduce the error until a tolerance is reached. 2. Write the error DOFVector in every adaption iteration to a file. 3. Calculate the exponent in the error estimate `\(\|u_h - I_h u^\ast\|_\#\leq C h^k\)` 4. (optional) Transform the mesh, by `\(x\mapsto 2x\)` and solve again. --- # Some hints ### Used functions/classes: ``` DOFVector << EXPRESSION; integrate( EXPRESSION ); // EXPRESSION: {valueOf(U), gradientOf(U), X(), X(i), +,-,*,/, // unary_dot(EXPR.), pow<2>(EXPR.), absolute(EXPR.)} RefinementManager* refManager = prob.getRefinementManager(); refManager->globalRefine(prob.getMesh(), NR_OF_REFINEMENTS); adaptInfo.getSpaceTolerance(0); adaptInfo.getMaxSpaceIteration(); io::writeFile( DOFVECTOR, FILENAME ); ``` ### Parameters to modify: ```matlab adapt[0]->tolerance: DOUBLE adapt->max iteration: INTEGER ```