Timeline for Solving Poisson equations as mixed Laplace using $RT_0-P_0$ pair
Current License: CC BY-SA 4.0
24 events
| when toggle format | what | by | license | comment | |
|---|---|---|---|---|---|
| Jan 8, 2022 at 23:30 | comment | added | bob_bill | Thanks a lot @knl, for the mathematical side I was already using those classics, while I have to say FEniCS book was new to me, but really cool (I already saw most of Bangerth's lectures notes). Your suggestions are really precious, thanks again :-) | |
| Jan 7, 2022 at 7:13 | comment | added | knl | I read so many books that I've lost count. There isn't a single source I used. However, reading books by people who have developed codes themselves and comparing them to some of the more mathematical texts have helped me a lot. Some examples are Demkowicz' books, Larson book is very nice although the codes are not very fast, Wolfgang Bangerth's lecture notes, FEniCS book is quite good overview of modern thinking. The mathematical side is: Braess, Johnson, Ern-Guermond books, Brenner-Scott, and Boffi-Brezzi-Fortin. However, I also learn a lot from colleagues in academia so hard to say. | |
| Jan 5, 2022 at 12:12 | vote | accept | bob_bill | ||
| Jan 5, 2022 at 12:12 | history | bounty ended | bob_bill | ||
| Jan 5, 2022 at 12:12 | comment | added | bob_bill | You're right. Thanks so much @knl, you've been really helpful. I accepted your answer. I've seen you're developing a python library. Is there any reference/book that describes FEM implementations from the Mathematical perspective (i.e. by thinking in terms of DoFs, not like engineers) that you've used to develop it, or is it "only experience"? | |
| Jan 5, 2022 at 7:34 | comment | added | knl | Note that even in standard $H^1$ FEM we are doing two transformations, the other one is just so trivial that you never think about it as a transformation. ($\psi(x) = \hat{\psi}(\hat{x})$) | |
| Jan 5, 2022 at 7:31 | comment | added | knl | Yes, we are applying two transformations. $F_K$ defines the shape of the global triangle. Piola makes sure that the derivative DOFs are preserved through the transformation. You can think about it in 1D. If you have $f(u) = u(a)$ as a DOF we can transform basis function as $\psi(x) = \hat{\psi}(\hat{x})$ and, as a consequence, $f(\hat{\psi}) = f(\psi)$. However, if you have $g(u) = u^\prime(a)$ as a DOF we would have $g(\hat{\psi}) \not= g(\psi)$. Hence, we cannot do $\psi(x) = \hat{\psi}(\hat{x})$ and need to invent something else. (Piola preserves normal derivatives in 2D and 3D). | |
| Jan 5, 2022 at 0:37 | comment | added | bob_bill | Sorry for asking again, but the "dx" term with $\det(B_K)$ (which simplifies with the $\det(B_K)$ in the divergence term) then will appear because we're still using the transformation $x = B_K \hat{x} + b_K$, i.e. $x=F_K(\hat{x})$ when we come back to $\hat{K}$, right? What is puzzling me is that is seems that we're applying two transformations: Piola and $F_K$. I understand the definition of Piola transformation, but I don't see how the change of variable from $K$ to $\hat{K}$ is working in this case. @knl | |
| Jan 4, 2022 at 21:24 | comment | added | knl | I'm not sure what is going on. The contravariant Piola transformation is defined as $$\frac{1}{\det(B_K)} B_K \hat{\Psi}(\hat{x}) = \Psi(x).$$ This implies that the divergence transforms as $$\frac{1}{\det(B_K)} \hat{\mathrm{div}}\hat{\Psi} = \mathrm{div}\Psi.$$ | |
| Jan 4, 2022 at 19:19 | comment | added | bob_bill | So basically we have $\int_K \operatorname{div}(\Phi)p(x) dx$ and since $x = F_K(\hat{x})$ and we $$\int_{\hat{K}} \operatorname{div}(\frac{1}{\det(B_K))} \hat{\Phi}(\hat{x})) \hat{q_j}(\hat{x}) |\det(B_K)| d \hat{x}$$ What is not clear to me is why we're allowed to put that $\frac{1}{\det(B_K)}$ inside the divergence, instead of the usual affine map. This is not mathematically correct, right? @knl | |
| Jan 4, 2022 at 14:48 | comment | added | knl | The code in your reference is a bit unreadable due to the unnecessarily complicated vectorization. However, it seems in the first snippet they multiply by detJ and in the second snippet they divide by detJ so the end result should be the same. The equation (9) is really confusing because the divergence has no "hat" implying that it would be in the global coordinates. Somehow I feel that shouldn't be the case so I assume they have a typo somewhere in the formulae. | |
| Jan 4, 2022 at 14:40 | comment | added | knl | Usually you multiply by detJ because of the integral transformation. The integral is transformed to the reference element and, as a consequence, we multiply by detJ. However, Piola transformation requires dividing by detJ so those two cancel out and that's why your approach works. | |
| Jan 4, 2022 at 10:58 | comment | added | bob_bill | You spotted the biggest typo. Thanks a lot. I edited my question with my results and the new code. It seems to work (see solution and EOC w.r.t DoFs), even if I do not understand what should be the local contribution to the matrix $B$, since in the paper I am using it's not clear to me what's going on. If you could clarify this last aspect, then I have no other questions (I opened a bounty so to give you more credits, as you looked also the code) | |
| Jan 2, 2022 at 22:09 | history | edited | knl | CC BY-SA 4.0 |
added 331 characters in body
|
| Jan 2, 2022 at 22:05 | comment | added | knl | Looking at your code, that could explain the error. What you need to do is to uniquely index all edges in your mesh and then use those instead of the vertex indices when assembling. I edited the answer to reflect this. | |
| Jan 2, 2022 at 22:05 | comment | added | knl | RT0 basis functions are one per edge, not one per vertex. The indexing is different. | |
| Jan 2, 2022 at 19:47 | comment | added | bob_bill |
Uhm, then there must be an error in my assembly routine. I have, in each triangle $K$, three DoFs for the velocity $\boldsymbol{u}$ whose global index is given by the l2g vector, and one for the pressure. On each velocity DoF (corresponding to the vertex of the triangle) I am assuming to have a RT basis function, so the "divergence" matrix is a 3 by 1 in each $K$, while the mass matrix is a 3by3. That's correct, right? @knl
|
|
| Jan 2, 2022 at 11:22 | comment | added | knl | No postprocessing unless you want to do Neumann boundary condition. | |
| Jan 2, 2022 at 10:26 | comment | added | bob_bill | So, correct me if I'm wrong, in my case I shouldn't post process the matrices $A,B$ after the assembly, right? The problem then must be in how I assemble $B$, since its nullspace is not made by a multiple of the vector $1$. Maybe the problem is in how I considered the DoFs for $B$ in the snippet @knl | |
| Jan 2, 2022 at 8:15 | comment | added | knl | Zero Dirichlet is a natural BC and automatically satisfied. | |
| Jan 1, 2022 at 16:37 | comment | added | bob_bill | Okay. For what concerns the homogeneous Dirichlet BC, how should I impose them? In my code, I just assembled the global block matrix and solved using that matrix, without any postprocess, because the term with the Dirichlet BC in the weak form is identically $0$. @knl | |
| Jan 1, 2022 at 15:27 | comment | added | knl | Alright, in the DOF you are dividing by the length of the edge. Then it's fine I guess. | |
| Jan 1, 2022 at 10:59 | comment | added | bob_bill | Honestly, I just found the expression for those functions and used them. I took them from siam.org/Portals/0/Publications/SIURO/Vol12/S01743.pdf (page 250) | |
| Jan 1, 2022 at 7:28 | history | answered | knl | CC BY-SA 4.0 |