5

Why would you want to do things on your own? The libraries you mention have all been run on 10,000+ cores and under the hood use PETSc, Trilinos, hypre, ... for the solution of linear systems or use matrix-free approaches. You would have to invest tens of man-years of work to implement the functionality and optimizations that has gone into these libraries -- ...


4

The reason that this particular mesh does not give the correct, uniform displacement solution to this problem is that it is "non-conforming." Specifically, at the intersection of the two cubes in the model, the two element edges that cross that face don't align with each other but instead cross each other. A typical face in a non-conforming mesh ...


4

You already know that at least theoretically, unconstrained matrices have a null space and consequently eigenvalues that are equal to zero. But, in practice, this is a meaningless condition because it can not be checked in an efficient way for large problems. The question you specifically ask is how you can detect whether constraints have been applied, and ...


3

The correct answer to your question, IMHO, is "depends on your target and your problem at hand". 1.) If your target is to simulate a large-scale problem on HPC and if you know of an existing code which can model the physics of your problem readily, then use the existing code. 2.) If an existing code does not yet support the physics of your problem ...


3

TL;DR How can I determine which constrain I need to apply to the system to make problem solved? Or how can I determine which rigid body constrain I should apply to the system? The constraints are given by the boundary conditions of your problem, so you should know them before you have a numerical method as the FEM. In that sense, that is more a physics ...


2

As you said, "If displacement not be constrained, equation above can not be solved, because the system can have rigid body motion" So you should try to apply constraints that will not allow the body to move i.e. translate or rotate. In 2D there are 2 translations (along x and y axis) and one rotation (along z axis) to be killed. In 3D there are 3 ...


2

The pressure is a discontinuous function, so you can't store nodal values at the vertices of the cells -- because then all neighboring cells would have the same value. Rather, if you want, you can pick any three points in the interior of the cell. In practice, people often choose three of the four vertices but consider these points logically part of the &...


2

The system of equations after discretization with the FEM can be written in real algebra as $$\begin{bmatrix} A_{R} &-A_{I}\\ A_{I} &A_{R} \end{bmatrix} \begin{Bmatrix}p_R\\ p_I\end{Bmatrix} = \begin{Bmatrix}f_R\\ f_I\end{Bmatrix}\, , $$ where the subindices $R$ and $I$ refer to the real and imaginary parts of the impedance matrix, pressure vector ...


1

I would highly suggest you go with an available FEM open-source library (say deal.II, FENICS, MFEM, etc.) instead of writing your own FEM code and then using PETSC as the underlying parallel algebra library. First, the majority of open source HPC FEM code already use either PETSC or Trilinos under the hood (deal.II supports both, FENICS uses PETSC, etc.). ...


1

As pointed out by @nicoguaro, I did some more tests with a scatterer. The overall set up of these new simulations is the same as in my previous answer, as well as the conclusion. There must be some way of proving mathematically that $A_I = 0$, but unfortunately I don't have the time to investigate it right now. I'll leave an image for further reference.


1

I'd like to share some tests I did in order to see if there's any penalty by not taking the imaginary part into account. Normally, the weak form of the Helmholtz equation is written as $$\int_{\Omega}(\nabla p \cdot \nabla \overline{q} - \kappa ^2 p \overline{q})d \Omega = \int_{\Gamma}g \overline{q}ds$$ But when we account for the real ($p_r$) and imaginary ...


1

To expand on Wolfgang Bangerth's answer, I think P0 DG schemes reduce to two-point cell-centered finite volume schemes. I don't know if DG convergence analysis always includes $p = 0$, but the resulting finite volume schemes can be shown to converge under appropriate "mesh orthogonality" conditions. https://math.unice.fr/~minjeaud/Donnees/...


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