Collapse vertex and apply constraint

[Daniel Garcia-Sanchez]

Hi,

I'm doing 3D electromagnetic simulations with Nelelec and NedelecSZ elements. I would like to collapse some vertices in other to follow the geometry.

The matrix should be singular if you collapse a vertex but according to Luca Heltai:

You can get away with it if you add enough constraints to your system so that each collapsed degree of freedom is identical (you can do this by using the ConstraintMatrix class). This should be enough to remove the singularity from the matrix (and it is often used as a trick to deal with singular integrals).

https://groups.google.com/forum/#!searchin/dealii/degenerate%7Csort:relevance/dealii/VLXWrUVh_18/ohiblJOlCtUJ

I gave it a try with a simple simulation (step-6) with 25 dofs and a rectangular mesh.

I used add_line() and add_entry() from AffineConstraints

If I collapse the vertex and don't apply the constraint, I still get the right solution and the determinant is non zero. I expected the determinant to be zero for this particular case. The matrix should be singular.

I did 3 test which gave the same solution:

• Original mesh. determinant(system_matrix) = 6.6e10
• Move vertex and do not apply constraint. determinant(system_matrix) = 6e11
• Move vert and apply constraint. det(system_matrix) = 1.2e12

Below you can find the simulations as you can see I always get the right solution for all the cases.

Why the determinant when I collapse the vertex and don't apply the constraint is non-zero? The matrix should be singular.

How do I find the dof asociated to a vertex? I just did try and error for this particular test.

Can I do the same trick in 3D for the Nedelec and NedelecSZ elements?

How do I check the quality of a mesh and that the matrix is not singular? I can't do the determinant for a large matrix.

Thanks!

Daniel

[Daniel Garcia-Sanchez]

Hi,

I realized that the mesh is not clear in the solution images. Below you can find the meshes: with and without a vertex collapse.

Thanks!

Daniel

[Wolfgang Bangerth]

On 7/30/19 1:28 PM, Daniel Garcia-Sanchez wrote:
>
> I did 3 test which gave the same solution:
>

> * Original mesh. determinant(system_matrix) = 6.6e10
> * Move vertex and do not apply constraint. determinant(system_matrix)
> = 6e11
> * Move vert and apply constraint. det(system_matrix) = 1.2e12

You can't compute the determinant for anything but quite small matrices.
It's not a numerically stable operation. The only way to determine
whether a matrix is singular is to compute all eigenvalues and see
whether one or more are "small" compared to the size of the other
eigenvalues.

Can you do that for the matrices you have, and plot the eigenvalues in
the same plot? Is there a qualitative difference?

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Wolfgang Bangerth]

[Daniel reports that his answers get blocked, so I'll post for him here, omitting pictures.]

Hi Wolfgang,

I calculated the eigenvalues and the condition number. It seems that
it is ok to collapse one vertex and transform the quad into a
triangle.

Interestingly, applying the constraints with AffineConstraint does not
change the condition number. Even if I fix all the dofs with
add_line() the condition number doesn't change.

I collapsed one vertex in a 3D simulation with NedelecSZ and the
solution and the condition number are ok. I have to do more tests with
NelelecSZ.

In the simple simulation from step-6. I collapse one vertex and then
if I move another vertex and transform the quad into a line the
surface becomes very small and the condition number becomes very
large. Note that I didn't completely squeeze the quad into a line. If
the surface is zero, then as expected the simulation breaks.

This means that it is ok to collapse a vertex as long as the
surface/volume does not become too small?
According to this publication it is ok to collapse vertices
https://www.sciencedirect.com/science/article/pii/S1877705814016713

Enclosed you can find my program and the python script to analyze the data.
To run it: ./collapse_vertex>matrix.txt

You can also find enclosed the meshes. I added the condition number in
the image for each mesh.

Below you can find the detailed results :

- Original grid:
* surface_ratio = 1
* eigenvalues = [  0.666667     0.666667     0.666667     0.666667
1.17200912   1.33333
   1.33333      1.33333      1.33333      1.33333      1.33333      1.33333
   1.33333      1.33333      1.33333      1.33333      1.33333      9.13881996
   9.13881996   9.42425     12.79450083  21.95213004  21.95213004  31.83336
  58.59408006]
* condition_number = 80

- Collapse vertex a
* surface_ratio = 0.5
* eigenvalues = [   0.666667      0.666667      0.666667      0.666667
     1.17122833
    1.33333       1.33333       1.33333       1.33333       1.33333
    1.33333       1.33333       1.33333       1.33333       1.33333
    1.33333       1.33333       8.92744864   10.04655013   12.5454928
   13.98137723   20.47958898   31.17725707   45.95448961  185.24586722]
* condition_number = 3e+02

- Collapse vertex a, move vertex b
* surface_ratio = 0.1
* eigenvalues = [   0.666667      0.666667      0.666667      0.666667
     1.16902972
    1.33333       1.33333       1.33333       1.33333       1.33333
    1.33333       1.33333       1.33333       1.33333       1.33939
    1.52754       1.5336        5.11586123   10.46272767   12.09280555
   14.54544349   22.42497592   31.57412136   50.19072911  264.78640594]
* condition_number = 4e+02

-Collapse vertex a, almost collapse vertex b
* surface_ratio = 1e-8
* eigenvalues = [  6.66667000e-01   6.66667000e-01   6.66667000e-01
6.66667000e-01
   1.21307648e+00   1.33333000e+00   1.33333000e+00   1.33333000e+00
   1.33333000e+00   1.33333000e+00   1.33333000e+00   1.33333000e+00
   1.33333000e+00   1.33333000e+00   1.35897000e+00   1.66666000e+00
   1.69231000e+00   8.83186491e+00   1.10689839e+01   1.37682911e+01
   1.98164659e+01   2.84090528e+01   4.79480606e+01   2.92923117e+05
   1.70718688e+06]
* condition_number = 3e+06

Best, 

Dani 

[Wolfgang Bangerth]

> I calculated the eigenvalues and the condition number. It seems that
> it is ok to collapse one vertex and transform the quad into a
> triangle.
>
> Interestingly, applying the constraints with AffineConstraint does not
> change the condition number. Even if I fix all the dofs with
> add_line() the condition number doesn't change.

Yes. Adding constraints is the same as adding information, which can only make
the condition number better. (Not a scientific answer, but intuitively true
anyway.)

> I collapsed one vertex in a 3D simulation with NedelecSZ and the
> solution and the condition number are ok. I have to do more tests with
> NelelecSZ.
>
> In the simple simulation from step-6. I collapse one vertex and then
> if I move another vertex and transform the quad into a line the
> surface becomes very small and the condition number becomes very
> large. Note that I didn't completely squeeze the quad into a line. If
> the surface is zero, then as expected the simulation breaks.

That's correct too -- if you just squeeze an element but not completely
collapse + add constraints, then that's equivalent to a mesh with a poorly
shaped cell. Since the condition number of the matrix involves the ratio of
the largest to the smallest eigenvalue of the Jacobian of the mapping on each
cell, you end up with a large condition number if you have poorly shaped
cells. This happens here:

> -Collapse vertex a, almost collapse vertex b
> * surface_ratio = 1e-8
> * eigenvalues = [ 6.66667000e-01 6.66667000e-01 6.66667000e-01
> 6.66667000e-01
> 1.21307648e+00 1.33333000e+00 1.33333000e+00 1.33333000e+00
> 1.33333000e+00 1.33333000e+00 1.33333000e+00 1.33333000e+00
> 1.33333000e+00 1.33333000e+00 1.35897000e+00 1.66666000e+00
> 1.69231000e+00 8.83186491e+00 1.10689839e+01 1.37682911e+01
> 1.98164659e+01 2.84090528e+01 4.79480606e+01 2.92923117e+05
> 1.70718688e+06]
> * condition_number = 3e+06

For this poorly shaped cell, you can compute the eigenvalues of the Jacobian
of the mapping, if you want. You'll find that they behave roughly like Delta x
and Delta y, respectively, and for your mesh, these are very different.