11

Julia is built in such a way that you will never see a full PETSc-like library, and that's on purpose. PETSc is not a single thing: it is an HPC library with some utility functions, linear solvers, nonlinear solvers, ... a whole can of soup that works within its own world but not outside of it. Julia's package ecosystem is build on generics and composability....


9

This paper by Kirby and Mitchell describes the implementation of $C^1$ elements in the Firedrake package*. One of the main use cases is biharmonic problems, which show up in the elastic deformation of thin plates, or other higher-order PDE like Cahn-Hilliard. My impression from reading the paper and talking to Rob Kirby is that $C^1$ elements are better but ...


9

You made an error in the indices for the real and imaginary part. $$ M\frac{\xi_{R,n+1}-\xi_{R,n}}{Δt}=-A\frac{\xi_{I,n+1}+\xi_{I,n}}2 \\ M\frac{\xi_{I,n+1}-\xi_{I,n}}{Δt}=A\frac{\xi_{R,n+1}+\xi_{R,n}}2 $$ has all 4 vectors in each equation, it does not factorize the way you wrote. $$ \begin{bmatrix} 2M&Δt\,A\\ -Δt\,A&2M \end{bmatrix} \begin{bmatrix} ...


9

Short answer: Just replicate the vector of interpolation functions into a block-diagonal matrix, as showed e.g. on page 5 in this lecture note. Detailed answer: Mathematically oriented texts typically do not bother about detailing the implementation, that is why you haven't met it. On the other hand, engineering FEM do detail the implementation (but lack the ...


8

$C^1$ elements are mostly a historic relic. In the finite element method, the traditional view is that the best methods are "conforming", i.e., methods where the finite element space $V_h$ is a subspace of the space $V$ in which the solution lies. For second-order elliptic equations, $V=H^1$ and functions that are continuous and piecewise ...


7

Yes. The constants that appear in the interpolation estimates upon which finite element error estimates are based contain minimum and maximum angles of triangles/tetrahedra (or similar geometric measures for quadrilaterals/hexahedra). These constants are smallest whenever you have equilateral triangles or square quadrilaterals.


7

Your reasoning for vector type equation is correct when you wrote it down for each individual component. What occurs is that you need to write your interpolation basis function and your test function as vectors. Hence, starting from: $$\nabla^2 \mathbf{u} =0 $$ in strong integral form is: $$\int \mathbf{v} \cdot \nabla^2 \mathbf{u} \, d\Omega = 0 $$ where $\...


7

The other answer has everything you already need, but it's also worth pointing out that $u_h$ is computable whereas $I_hu$ is not: The latter requires you to know the exact solution, which in general we of course don't know (and if we did, we didn't need to compute a Galerkin approximation $u_h$).


6

In my opinion the basics stay the same but there are some differences worth commenting upon. The basic principles of defining a local coordinate system, integrating local dense contribution matrices using quadrature rules, and assembling them into global sparse matrices are identical. Often, hexahedral basis functions arise from taking tensor products of ...


6

Let's assume $u$ is the solution to the variational problem and $u_h$ is the Galerkin approximation of the solution on the subspace $V_h \subset V$. By Cea's lemma you have: \begin{equation} ||u-u_h||_V \leq C\inf\limits_{v_h \in V_h}||u-v_h||_V \end{equation} for some positive constant $C$. Now you define the projection $I_h:V \rightarrow V_h$. As you ...


6

I am biased, because I'm one of the authors, but would suggest looking at the deal.II library. It has an extensive tutorial.


6

I realized I made a mistake in the question by claiming that the reverse multiplication order does not yield the same result. The correct answer is as follows. My main quest is to find the following components (needed for calculating the second derivatives of shape functions): \begin{equation} \frac{\partial r}{\partial x}, \quad \frac{\partial r}{\partial y}...


6

You really don't want to implement this yourself -- you'll spend a year or two on things others have already done, and will have done far better than you can hope for. The difficulty is generally getting h and p refinement to work at the same time. That is not a trivial challenge. The implementation in deal.II is largely described here: https://www.math....


5

As a preamble, I would not expect that splitting $E$ into real/imaginary parts is very profitable. Normally, block 2x2 systems are motivated because one block of unknowns is "easier" to solve than the other in some sense (better conditioned? smaller in cardinality? etc). This is not the case for time-harmonic Maxwell, I think you'd be better off ...


5

It turns out that I have just the right paper for you on this subject: https://www.math.colostate.edu/~bangerth/publications/2013-pattern.pdf


5

The triangle inequality is your friend. Let's ignore the issue of boundary approximation for a moment, then you are computing a solution with inexact linear solver. Let's call it $u_h$. We will call the exact solution of the PDE $u$, and the exact finite element solution $u_\text{FEM}$; neither of these can be computed exactly (in the case of $u_\text{FEM}$ ...


5

$\vec\Phi = K_i \nabla u_i$ is the flux across the interface. For example, if $u$ is the thermal energy density and $K$ the thermal conductivity, then $\vec\Phi$ is the thermal energy flux. Energy conservation then dictates that whatever flows into the interface on one side ($\vec n_1 \cdot K_1 \nabla \vec\Phi_1$) better be equal to what flows out on the ...


4

$h$ is a measure of the mesh size. In the example, they are using rectangular elements. For which a commonly used measure of the mesh size is the length of the largest diagonal. Looking at the table above the sentence you quoted: $(M_5-M_6)/(M_4-M_5)=(0.14052586 −0.14056422)/(0.14037251−0.14052586)\approx 1/4 = (h_6/h_5)^2$ where $M_i$ is the mean at $i$-th ...


4

For both schemes: I am going to lump $B$ and $C$ into a single matrix $B$ for convenience of communicating the idea. Likewise, $y$ and $z$ are now $y$, and $F_y$ and $F_z$ are $F_y$. $Ax+B^Ty=F_x$ or $x = A^{-1}(F_x-B^Ty)$ $Bx=F_y$ or $BA^{-1}(F_x-B^Ty)=F_y$ $BA^{-1}B^Ty=BA^{-1}F_x-F_y$ Note that $BA^{-1}B^T$ is symmetric and invertible, and can be ...


4

You normally don't assign material properties to a node. They are part of a subdomain and, hence, "assigned" to an element. They appear in the integrals you compute, e.g. $$\int\limits_{\Omega_e} k \nabla u \nabla w d\Omega\, ,$$ here $\Omega_e$ represents the element identified as $e$. Nevertheless, due to interpolation the computed integrals are ...


4

One approach to convert this into an ODE is with index reduction methods. These allow you to convert high-index DAEs into low-index DAEs or ODEs. See section VII.2 of "Solving Ordinary Differential Equations II" by Hairer and Wanner. Consider a generic, Hessenberg index-2 DAE $$ \begin{align} y' &= f(y, z) \\ 0 &= g(y) \end{align} $$ ...


4

You can follow the concept of Pascal's pyramid (see the picture) for identifying the basis functions for elements of different orders. The picture is taken from what-when-how.com. Refer to the linked site for the details on computing the coefficients.


4

Let's take an $n\times n$ mesh (with $N=n^2$ unknowns) and think about whether you can enumerate them in such a way that you end up with a bandwidth less than $m=n=\sqrt{N}$? You get this bandwidth with a 5-point stencil if you enumerate the first row left to right, then the next row left to right, etc. In that case, each degree of freedom $i$ couples with $...


4

The following paper is a good starting point. Arnold, D. N. (1990). Mixed finite element methods for elliptic problems. Computer methods in applied mechanics and engineering, 82(1-3), 281-300. Author copy: https://www-users.cse.umn.edu/~arnold/papers/mixed.pdf It presents the general problem and why one would like to use one formulation or the other. Also, ...


4

Why not try drawing the basis? You use Raviart-Thomas basis to approximate a vector field so a quiver plot makes sense. This basis function (one out of three) is associated with the left edge. It is a kind of flux towards the left edge which is orthogonal to the normals of the other edges and magnitude is zero on the opposite vertex. In the global basis you ...


4

In addition to Wolfgang Bangerth's explanation of temperature and concentration, let me give an other application where such interface conditions arise: (linear) elasticity, which has a similar structure to the elliptic equation in your example. Consider the situation with no interface first. Then the equilibrium, constitutive and kinematic equations of ...


4

From an old set of notes, I think that since your material is incompressible, the stress is determined by the strain energy density function $W$ only upto the hydrostatic pressure. The Cauchy stress is given by $$ \pmb{\sigma} = -p\pmb{1} + 2\pmb{F}\frac{\partial{W}}{\partial{\pmb{C}}}\pmb{F}^T $$ The second Piola-Kirchhoff stress in incompressible ...


4

We can show how it works on the example of linear elasticity. In classical finite elements formulations, on every node, we will have a scalar shape (base) function to which we have associated number coefficients equal to the dimension of the problem. Coefficients at nodes are interpreted as physical nodal displacements. Let displacements are are approximated ...


3

Here is at least an idea, whether it works is a different question. Let's say you sort unknowns so that you have the ones in the interior of the domain first, and then all those at the boundary. Then the matrix that corresponds to your problem decomposes in the following way: $$ A = \begin{pmatrix} A^{\circ\circ} & B^{\circ\partial} \\ C^{\...


3

Using the typical expansion functions (1-forms/edge-elements for E, and 2-forms/facet-elements for B) the formulations are basically the same after spatial discretization and you'd expect more or less the same accuracy. I do think they express slightly different opinions about time integration. The mixed E/B formulation nudges you in the direction of ...


Only top voted, non community-wiki answers of a minimum length are eligible