0
$\begingroup$

I asked this question on the math stack exchange and got an answer, but I am just as utterly confused as before. My fundamental goal is to actually construct the matrix, that is, a series of steps I can tell an actual computer to build it up.

The gist of it is that there is a 2004 paper called "An Intuitive Framework for Real-Time Freeform Modeling". I am interested in equation (5) on the paper.

I am trying to build that matrix, but it is not obvious to me from the paper alone what the actual coefficients need to be. The answer that I got in the linked post seems correct but it spends most of the time talking about the derivation of the specific matrix, but right now I don't care about the FEM theory used to derive the matrix, all I am trying to do is implement the algorithm. I just want to know what coefficients go where.

I know that many of the rows of $L$ (matrix in the paper) will have the coefficients specified by the discrete Laplace-Beltrami operator.

My current understanding is as follows:

Assume that your vertices are sorted such that $k = m + n$ where $k$ is the total number of vertices, $m$ the number of free vertices and $n$ the number of fixed vertices. We will assume that $i < m$ means that the vertex $v_i$ is free, and fixed otherwise.

The way to construct the matrix $L$ for a single level (i.e. the left most image in figure 2) would be:

for vertex in free_vertices:
    for neighbour in beighbours of vertex:
        wij = cotan_weights(vertex, neighbour)
        W[vertex, neighbour] = wij 
        W[vertex, vertex] += -wij 
    area = vertex_area(vertex)
    M[vertex, vertex] = 1.0 / vertex_area

for vertex in fixed_vertices:
    M[vertex, vertex] = 1.0 
    W[vertex, vertex] = 1.0

And then the full Laplacian is $L = M \times W$, and it has dimensions $k\times k$. Is this at all correct? If this is incorrect, what is the right way to construct the matrix?

Additional info

I am trying to get this to work in the simplest possible example I could think of, this simple small cylinder:

enter image description here

In this case the top and bottom vertices are fixed the middle are allowed to move.

In this case I get the following matrix, if I try it with my formulation:

-4.2, 0.7, 0.0, 0.7, 0.0, 0.0, 0.0, 1.4, 0.0, 0.0, 0.0, 1.4, 
0.7, -4.2, 0.7, 0.0, 0.0, 0.0, 1.4, 0.0, 0.0, 0.0, 1.4, 0.0, 
0.0, 0.7, -4.2, 0.7, 0.0, 1.4, 0.0, 0.0, 0.0, 1.4, 0.0, 0.0, 
0.7, 0.0, 0.7, -4.2, 1.4, 0.0, 0.0, 0.0, 1.4, 0.0, 0.0, 0.0, 
0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0,

Which is yielding incorrect results:


x in: [0.0, 0.0, 0.0, 0.0, 0.0, -1.0, -0.0, 1.0, 0.0, -1.0, -0.0, 1.0]
x out: [0.0, -21.0, 21.0, 0.0, 0.0, -1.0, -0.0, 1.0, 0.0, -1.0, -0.0, 1.0]
y in: [0.0, 0.0, 0.0, 0.0, -1.0, -0.0, 1.0, 0.0, -1.0, -0.0, 1.0, 0.0]
y out: [28.0, -44.2, 3.5, -0.7, -1.0, -0.0, 1.0, 0.0, -1.0, -0.0, 1.0, 0.0]
z in: [0.0, 0.0, 0.0, 0.0, 2.0, 2.0, 2.0, 2.0, 0.0, 0.0, 0.0, 0.0]
z out: [-28.0, 23.3, -24.5, 0.7, 2.0, 2.0, 2.0, 2.0, 0.0, 0.0, 0.0, 0.0]

I know that libigl implements this in example 401]4, which is exactly what i am trying to achieve. But unfortunately the implementation of the harmonic subroutine is very impregnable and I have not been able to reverse engineer what the matrix is actually supposed to look like.

Improve this question
$\endgroup$
2
  • $\begingroup$ How did you decide that your solution is incorrect, because of the large values? Note that your problem is nonlinear. As far as the definition of harmonic goes, it does: $Q = -L$, and then has a loop that concatenates $Q = Q M^{-1}(-L)$. So you basically get $Q = -L(-M^{-1}L)^{k-1}$, $L$ being the cotangent stiffness matrix and $M$ being the mass matrix. Then instead of specifying the boundary conditions explicitly as you do, they use a Lagrange multipliers approach derived from the problem: $\min_u u^TQu, \,:\, Cu = Cb$, which yields $\max_{\lambda}\min_u u^TQu +\lambda (Cb-Cu)$. $\endgroup$
    – lightxbulb
    Commented Sep 21, 2023 at 11:34
  • $\begingroup$ The matrix $A$ is a function of position, so you actually have a non-linear problem. To solve it you can instead try to find the steady state of the time-evolution equation $\partial_t p = \Delta_S p$ with Dirichlet BCs. Discretising time and space in a semi-implicit manner yields: $\frac{M^kp^{k+1}−M^kp^k}{\tau}=L^kp^{k+1}$, thus you can now consecutively solve the problems $(M^k−\tau L^k)p^{k+1}=M^kp^k$ (for the bc indices keep $p^{k+1} = p^k$). You can try a stopping criterion $\|p^{k+1}−p^k\|<\epsilon$. Note that the superscripts denote iterate indices and not powers. $\endgroup$
    – lightxbulb
    Commented Sep 21, 2023 at 13:18

1 Answer 1

Reset to default
3
+50
$\begingroup$

I checked your solution, it's wrong. I ran the following MATLAB code:

A = zeros(12,12);
A(1,:) = [-4.2, 0.7, 0.0, 0.7, 0.0, 0.0, 0.0, 1.4, 0.0, 0.0, 0.0, 1.4];
A(2,:) = [0.7, -4.2, 0.7, 0.0, 0.0, 0.0, 1.4, 0.0, 0.0, 0.0, 1.4, 0.0];
A(3,:) = [0.0, 0.7, -4.2, 0.7, 0.0, 1.4, 0.0, 0.0, 0.0, 1.4, 0.0, 0.0];
A(4,:) = [0.7, 0.0, 0.7, -4.2, 1.4, 0.0, 0.0, 0.0, 1.4, 0.0, 0.0, 0.0];
for i=5:12
    A(i,i) = 1;
end

x_in = [0.0, 0.0, 0.0, 0.0, 0.0, -1.0, -0.0, 1.0, 0.0, -1.0, -0.0, 1.0]';
y_in = [0.0, 0.0, 0.0, 0.0, -1.0, -0.0, 1.0, 0.0, -1.0, -0.0, 1.0, 0.0]';
z_in = [0.0, 0.0, 0.0, 0.0, 2.0, 2.0, 2.0, 2.0, 0.0, 0.0, 0.0, 0.0]';
x_out = (A\x_in);
y_out = (A\y_in);
z_out = (A\z_in);

The resulting solution is:

x_out = [0.666666666666667, 0,  -0.666666666666667, 0,  0,  -1, 0,  1,  0,  -1, 0,  1]
y_out = [0, 0.666666666666667,  0,  -0.666666666666667, -1, 0,  1,  0,  -1, 0,  1,  0]
z_out = [1.00000000000000,  1.00000000000000,   1.00000000000000,   1.00000000000000,   2,  2,  2,  2,  0,  0,  0,  0]

This clearly doesn't match what you have. It is very likely that your solver has a bug or doesn't work for the matrix above.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.