24
When I studied the finite element method in graduate school, this notion of multiplying by a weight function was also very alien to me. Eventually, I did find a nice (albeit non-rigorous) analogy that helped me understand it. This analogy is based on 3D vector geometry and an understanding of projections and dot-product.
3D Geometry
Imagine a 2D plane ...
5
The problem that you are seeing is a well known problem of these basis functions:
$$\phi_{m,n}(x,y) = x^{m}y^{n}$$
I quote from here:
The reason why the coefficient matrix is nearly singular and ill-conditioned is that our basis functions $\Psi_{i}(x) = x^{i}$ are nearly linearly dependent for large i. That is, $x^{i}$ and $x^{i+1}$ are very close for i ...
4
The projection scheme is indeed analytic, provided that the initial conditions are incompressible. You can derive it by taking the divergence of the evolution equation for $u$.
Finding the projection P can be tricky; it's easier with some methods like spectral / pseudospectral methods. But one has to be careful that numerical error doesn't move the ...
4
If you use a sharper trace inequality you get
$$\begin{aligned}\|u-Pu\|_{0,\partial \Omega}^2 &\lesssim \|u-Pu\|_0 \|u-Pu\|_1 \\
&\lesssim h^{k+1} h^{k}|u|^2_{k+1} \end{aligned} $$
implying
$$\|u-Pu\|_{0,\partial \Omega} \lesssim h^{k+\frac12},$$
just as Wolfgang suggested in the comments.
4
What about simply computing the projection of $w$ along the complement of $V$?
$$
\tilde w = [I - V(VV^T)^{-1}V^T]w,
$$
where $V$ is the matrix that has the basis vectors of $V$ as columns. This $\tilde w$ has all the properties that you look for.
If you are in the non standard scalar product, then there is a symmetric positive matrix $M$ that induces this ...
4
It's important to realize that the original velocity field $\mathbf v_0$ is also not divergence free in a pointwise sense. Rather, it is only divergence free when tested with the pressure test functions, i.e., $(\nabla\cdot\mathbf v_0,q)=0$ for all $q\in P_1(T_0)$ where $T_0$ is the original mesh.
What you want to achieve is to get a velocity field $\mathbf ...
3
In order to understand these results, you need to know how minimization and projection problems are connected. Namely,
Let $\mathbb V$ be a subspace of $\mathbb C^n$ and take $y \in \mathbb C^n$; then $\hat y = \text{argmin}_{x \in \mathbb V} || y - x ||$ iff $(\hat y, x) = (y, x)$ for all $x \in \mathbb V$.
Here $(.,.)$ is a scalar product (not ...
3
See this recent paper on an extension of stochastic gradient descent that could be used on your problem:
https://arxiv.org/abs/1511.03760
You could also apply Dykstra's algorithm (or any other algorithm that does alternating projections on convex sets) by setting a target value for $\mu^{T}c \leq \gamma $ and reducing it once feasibility has been achieved.
3
Here are some elements of answers to the three questions and references to alternative methods for spherical parameterization:
1. How to compute a voxelization of a given model ?
What it means:
It means embedding your surface into a 3D voxel grids, and determining which voxels have an intersection with some triangles of the surface.
How to do that:
...
3
Note: I originally posted an answer that I was not 100% pleased with, so I have revised it heavily. 1/20/2017
The projection method is not in an exact approximation to the full system in general. There are a few reasons for this, but the most immediate has to do with the boundary conditions. Let's look at the continuous setting first. Consider the unsteady ...
3
Let's say you want to solve the Laplace equation, $-\Delta u = f$. Ideally, of course, you'd like to find a function $u$ so that the residual is zero: $r(u) = -\Delta u - f = 0$. But $u$ is an infinite dimensional object which in general we cannot represent on computers, so we have to find finite dimensional approximations $u_h$. Since $u_h$ is not the exact ...
3
The answer to your first question is no -- the active set does not only grow, it can also shrink, even if your objective function is convex. An example is if you start at some point on the boundary, then the first step moves along one active constraint to find the minimum of the objective function along this segment; but then, the second step will likely ...
2
The projection you consider is local on every cell. With the spaces you cite, your condition 1 requires test functions of degree $k-1$ on each edge (constant, linear, and quadrature along the edge); for $k=3$, there are 3 linearly independent such shape functions per edge, for a total of 9. Your condition 2 requires test functions of degree $k-2$ for the ...
2
While this question is old and has been answered by plenty of smart people, I just want to jot down the intuition I use to explain the Galerkin method to people.
The goal in our situation is to find as close of a solution as we can to some continuous residual equation:
$$r(u) = 0$$.
Let us define the $i^{th}$ basis function as $\phi_{i}(x)$, define the ...
2
Boris Grigoryevich wants you not to be able to create residuals with the same functions you used to create the solution.
2
What you are observing is odd-even decoupling between the velocity and pressure. This is a well-known problem on collocated grids. At any point on the grid, the gradient of the pressure is calculated using $p_{i+1}$ and $p_{i-1}$:
$$\frac{u_i^{n+1}-u_i^n}{\Delta t} + (u_i^n.\nabla)u_i^n = -\frac{p_{i+1}-p_{i-1}}{2\Delta x}+\nu\nabla^2 u_i^n$$
Hence the ...
2
Converted to an answer from my comments to Jan's answer.
To fix dimensions and notation, let us say that $V$ is a $m\times n$ matrix with the vectors $v_i$ as columns.
As Jan notes, we have to compute
$$
\tilde w = [I - V(V^TMV)^{-1}V^TM]w,
$$
where $M$ is the symmetric $m\times m$ matrix associated to the scalar product.
If it were possible to solve this ...
2
In the matrix case, the least square solution of an overdetermined inconsistent system $Ax=b$ requires the solution of $(A^TA)x=A^Tb$. The projection $A(A^TA)^{-1}Ab$ reduces $b$ to the closest vector that is a consistent right side.
As in the matrix case, the solution of the overdetermined inconsistent system $∇z=V$ is reduced to a solvable problem by ...
2
Some thoughts on a particular case of low-rank matrices ($k=\text{rank}(M)\ll n)$. Here, I can suggest some economical version of SVDs:
Rand-SVD
Here is the R package documentation that also contains references for canonical papers of V. Rokhlin and P. Martinsson on this topic. That should reduce the complexity to $\mathcal O(nk^2+k^3)$ as opposed to $\...
2
Fortunately, LAPACK provides routines to deal with the $\mathbf Q$ factor from the $\mathbf A = \mathbf Q \mathbf R$ decomposition, [dgeqrf]. To find the projection of an arbitrary $\mathbf B$ onto the space orthogonal to $ \mathrm {range}(\mathbf A)$, you want to form $\mathbf C = \left(\mathbf I - \mathbf Q \mathbf Q^T\right) \mathbf B$. Here are two ...
2
The smooth solution turned out to have BC's applied in the following way:
Walls and inlet: $\frac{\partial p}{\partial n}=0$
Outlet: $p=0$
1
Thinking about it even more, I feel like the problem resembles triangulation of lines under orthographic camera model. It is possible and useful to visualize the setting in 3D:
The blue vector is what we like to recover. The identity element for the orthographic camera and projection matrices are defined respectively as:
$$
K=\begin{bmatrix}
1 & 0 & ...
1
The problem seems quite ill posed so I will assume some level of over-determinism is there.
First, note that your vector, $\mathbf{z}$ lives on a hypersphere. Thus, it is indeed on a $(k-1)$ dimensional manifold, $\mathcal{S}^{k-1}$. Then, because you don't know the scales, we could always write the problem as linear combination of $k$ unit vectors, $\{\...
1
Section 2.1 in the paper explains that PCA takes a list of vectors
$X = [X^{(0)} | X^{(1)} | ... | X^{(N-1)}]$
and maps it to
$T = [T^{(0)} | T^{(1)} | ... | T^{(N-1)}]$
So I presume $T$ contains the data in the new coordinate system.
It then explains that
$X = T P^T$
So I gather that $P$ contains the principle components.
I'm using the published ...
1
What sort of field is it? You've said "divergence free", but do you mean harmonic (zero divergence and curl), solenoidal (non-curl), or a mix of both? The distinction is important because it will affect the complexity of the algorithm. The goal is to do a Helmholtz decomposition numerically, interpolate the potentials, and then recalculate the velocities ...
1
Presuming you have a suitable solver (CPLEX, GUROBI, MOSEK, SCIP, many others), you can solve this as a Mixed-Integer Quadratic Program (by squaring the objective) or as a Mixed-Integer Second Order Cone Problem.
Here is code to solve it in YALMIP under MATLAB.
v = intvar(size(v_o))
optimize([-k <= v <= k, sum(v) == l],norm(v-v_0))
disp(value(v)) $ ...
1
tldr: Reformulate the projection and avoid the need for boundary conditions on the pressure.
I think you are misinterpreting the projection scheme. In all formulations that I know, the pressure is never really computed.
It rather goes like:
Compute a tentative velocity $u^*$ approximating $u^{n+1}$
Project $u^*$ onto the space of divergence-free functions ...
1
The most obvious choice is setting $\mathbf{u}^* = \mathbf{u}^{n+1}$ on the boundary.
This is true if your velocity is described by Dirichlet BCs, however, if your velocity BCs are described by Neumann BCs (fully developed flow, e.g.), you should use $\frac{\partial \mathbf{u}^*}{\partial n} = 0$. This, it sounds to me, is what you're looking for.
1
In the interior of cells, the Raviart-Thomas functions are continuous. As a consequence, the normal component is of course also continuous and the jump is zero. That may not be the case at places where the interface intersects cell boundaries, though. The set of these intersection points, however, consists of a finite number of points along a curve in 2d, or ...
1
That you get really big numbers is no accident. Consider
$$
C = \left(\begin{array}{ccc}
2 & 1 & \epsilon\\
1 & 2 & 1\\
\epsilon & 1 & 2
\end{array}\right)
$$
and
$$
X = \left(\begin{array}{ccc}
-1 & -1 & 1\\
-1 & 1 & 1\\
1 & 1 & 1\\
\end{array}\right).
$$
Note that $X$ has eigenvalues $\{1, \pm 2\}$ and $C$ ...
Only top voted, non community-wiki answers of a minimum length are eligible