Hot answers tagged

8

As one of the library's authors, I would of course love for deal.II to come out on top with this comparison. But I suspect it may not, and the answer lies in a factor you omit from your comparison: how long it actually takes to implement your code. Few people in academia with the skills to implement a FEM code from scratch spend more time solving PDEs than ...


5

Having written this entry, let me also answer you question :-) The issue in essence boils down to the following: given a function $u(x)$ and its interpolation $u_h(x)$ onto either a triangular or quadrilateral mesh with the same number of unknowns, then which of the two is more accurate? In other words, is $$ \| u - I_h^\text{triangles} u \| < \| u - ...


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

The SEM is a FEM! It's almost like all these different names are designed to confuse the newcomer. I will speak primarily about the most popular form which uses a tensor product Lagrange basis with nodes at the Legendre-Gauss-Lobatto (LGL) points. Once you have chosen that basis and integrate with the LGL quadrature using the same nodes. You have a method ...


4

The main advantage is that it reduces the Runge phenomenon and leads to faster convergence rates. It also presents less numerical dispersion and need less nodes per wavelength (see 1 and 2). So, I would say that you would prefer the method for wave propagation scenarios. Regarding software that includes SEM, I am aware of the following: FSELib: Matlab ...


4

Specific answers to this question are probably time-limited. However, the following general approach (from the great Eric S. Raymond) works very well: Rule of Modularity: Write simple parts connected by clean interfaces. Rule of Clarity: Clarity is better than cleverness. Rule of Composition: Design programs to be connected to other programs. Rule of ...


4

Unfortunately there's no tool for this. You can run each on a variety of input sizes to establish the computational complexity they appear to have, i.e. the $f$ in the $O(f(n))$ that characterizes each code. This can point you towards what underlying algorithms each is using and verify or not that something that should be $O(n)$ is actually implemented to ...


3

I fully agree with Wolfgang's answer, but I'd like to add something of my experience. I had recently to choose between tri and quad based FEM and after googling for it, I decided that I would code both from scratch for a fair comparison (I had a previous experience in writing FEM kernel. I used C++ and the Eigen library for matrices and vectors. I considered ...


2

I will try answer based on my experience with deal.ii. The max_per_row=5 means that at most there will be 5 non-zeros per row in the matrix. Since we now know this then we do not need to have a $1\times{}100$ matrix but rather a $1\times{}50$. In other words this parameter sets an upper-bound on the memory needed. In reality it is not stored as one vector ...


2

You go the grouping in the term $\nabla u^+-\nabla u^-$ mostly right but not quite. To see this, remember that the functions $\phi_i^{\pm}$ are discontinuous, and that consequently $u_0^+$ is not necessarily equal to $u_3^-$ -- if they were, you'd end up with a function that is continuous at that point. Instead, you need one degree of freedom per adjacent ...


1

You need to compute the projection $\Pi(f-cU_h)$ as a first step. This is a finite element field, so you can create a global field (in deal.II lingo: Create the corresponding finite element and a DoFHandler object) and compute it by global projection by inverting a mass matrix. Alternatively, because the field is discontinuous, you can also forgo the global ...


1

As mentioned by @Wolfgang Bangerth , the value may be at the interior points. However, to visualise that better in paraview, You can use a filter called "plot over line" which plots the data points over a given line, which could come very handy for these kinds of clarifications in future.


1

@cfdlab's answer gives a general way to construct solutions for discontinuous coefficients, but here is a slightly more theoretical perspective to it as well. All of this can be understood in 1d, so imagine all of the functions below to be functions of just one argument $x$ for a moment. First, you can't just choose any $u$ and $\alpha$ to obtain a function $...


1

If you want to use the coefficient of the form $$ \alpha = \begin{cases} \alpha_1 & r < 0.5 \\ \alpha_2 & r > 0.5 \end{cases} $$ it is useful to work in polar coordinates. You build two pieces of the solution $$ u_1(r,\theta), \qquad r < 0.5 $$ and $$ u_2(r,\theta), \qquad r > 0.5 $$ Then at the interface you ensure solution and flux ...


1

This is not really a program as much as a big repository of codes. John Burkardt from Florida State University maintains a rich collection of scripts in C++, Matlab and Fortran for a large range of problems. You can check out the FEM 2D library or go up a level and look at other codes listed there. These have good commenting on them and you can look at the ...


1

If a homogeneous Neumann problem is considered, i.e. $f=0$ in $\Omega$ and $\mathbf{n}\cdot\nabla\phi$ on $\Gamma$, one solution is given by the constant function $$ \phi\left(\mathbf{x}\right)=\phi_0\,. $$ Then the boundary integral equation reads $$\int\limits_{\Omega}\phi\left(\mathbf{x}' \right)\delta\left(\mathbf{x}'-\mathbf{x}\right)\mathrm{d}V' =\...


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