On Tuesday, July 30, 2019 at 10:06:59 PM UTC+2, Wolfgang Bangerth wrote:
>
>
> 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?
>
>
>
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 AffineConstraints does not
change the condition number. Even, if I fix all the dofs with add_line()
the condition number does not change.
I collapsed one vertex in a 3D simulation with NedelecSZ and the solution
and the condition are ok. I have to do more tests with 3D simulations with
NelelecSZ.
In the simple simulation from step-6. I collapsed one vertex and then, if I
move another vertex and squeeze the quad into almost a line the surface
becomes very small and then 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 in quads/hexes
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
[image: original_grid.png][image: collapse_a.png][image:
collapse_a_move_b.png][image: collapse_a_close_to_collapse_b.png]
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/c680299c-8124-4b8a-94fd-614a6e8f0017%40googlegroups.com.
#!/usr/bin/env python3
import numpy as np
import numpy.linalg
matrix_size = 25
matrix = np.zeros((matrix_size,matrix_size))
hola_txt = open('build/matrix.txt','r')
for line in hola_txt:
if line[0] == '(':
value = float(line.split(' ')[1])
coordinates = line.split(' ')[0]
x = int(coordinates.split(',')[0][1:])
y = int(coordinates.split(',')[1][:-1])
matrix[x,y] = value
determinant = np.linalg.det(matrix)
eigenvalues, eigenvectors = np.linalg.eig(matrix)
condition_number = np.linalg.cond(matrix)
print('{:.2e}'.format(determinant))
print(np.sort(eigenvalues))
print('{:.2e}'.format(condition_number))
/* ---------------------------------------------------------------------
*
* Copyright (C) 2000 - 2019 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE.md at
* the top level directory of deal.II.
*
* ---------------------------------------------------------------------
*
* Author: Wolfgang Bangerth, University of Heidelberg, 2000
*/
#include <deal.II/base/function_lib.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/vector.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <fstream>
#include <deal.II/fe/fe_q.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/numerics/error_estimator.h>
using namespace dealii;
template <int dim> class Step6 {
public:
Step6();
void run();
private:
void setup_system();
void assemble_system();
void solve();
void refine_grid();
void output_results(const unsigned int cycle) const;
Triangulation<dim> triangulation;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
AffineConstraints<double> constraints;
SparseMatrix<double> system_matrix;
SparsityPattern sparsity_pattern;
Vector<double> solution;
Vector<double> system_rhs;
const bool move_vertex = true;
const bool add_constraint = false;
const double displacement_vertex_a = 0.5;
const double displacement_vertex_b = 0.499999;
};
template <int dim> double coefficient(const Point<dim> &p) {
if (p.square() < 0.5 * 0.5)
return 20;
else
return 1;
}
template <int dim> Step6<dim>::Step6() : fe(1), dof_handler(triangulation) {}
template <int dim> void Step6<dim>::setup_system() {
dof_handler.distribute_dofs(fe);
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
VectorTools::interpolate_boundary_values(
dof_handler, 0, Functions::ZeroFunction<dim>(), constraints);
if (add_constraint) {
// number 5 is 8 counting from below
// number 8 is the center
unsigned int line_number = 8;
unsigned int other_dof = 5;
constraints.add_line(line_number);
constraints.add_entry(line_number, other_dof, 1);
}
if (add_constraint) {
// 13 is top
// 10 is below
// done 9,11,12,13
unsigned int line_number = 10;
unsigned int other_dof = 13;
constraints.add_line(line_number);
constraints.add_entry(line_number, other_dof, 1);
}
constraints.close();
DynamicSparsityPattern dsp(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints,
/*keep_constrained_dofs = */ false);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
}
template <int dim> void Step6<dim>::assemble_system() {
const QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators()) {
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit(cell);
for (unsigned int q_index = 0; q_index < n_q_points; ++q_index) {
const double current_coefficient =
coefficient<dim>(fe_values.quadrature_point(q_index));
for (unsigned int i = 0; i < dofs_per_cell; ++i) {
for (unsigned int j = 0; j < dofs_per_cell; ++j)
cell_matrix(i, j) +=
(current_coefficient * // a(x_q)
fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
fe_values.JxW(q_index)); // dx
cell_rhs(i) += (1.0 * // f(x)
fe_values.shape_value(i, q_index) * // phi_i(x_q)
fe_values.JxW(q_index)); // dx
}
}
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(
cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
}
}
template <int dim> void Step6<dim>::solve() {
SolverControl solver_control(1000, 1e-12);
SolverCG<> solver(solver_control);
PreconditionSSOR<> preconditioner;
preconditioner.initialize(system_matrix, 1.2);
solver.solve(system_matrix, solution, system_rhs, preconditioner);
constraints.distribute(solution);
}
template <int dim> void Step6<dim>::refine_grid() {
Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
KellyErrorEstimator<dim>::estimate(
dof_handler, QGauss<dim - 1>(fe.degree + 1),
std::map<types::boundary_id, const Function<dim> *>(), solution,
estimated_error_per_cell);
GridRefinement::refine_and_coarsen_fixed_number(
triangulation, estimated_error_per_cell, 0.3, 0.03);
triangulation.execute_coarsening_and_refinement();
}
template <int dim>
void Step6<dim>::output_results(const unsigned int cycle) const {
{
GridOut grid_out;
std::ofstream output("grid-" + std::to_string(cycle) + ".eps");
grid_out.write_eps(triangulation, output);
}
{
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "solution");
data_out.build_patches();
std::ofstream output("solution-" + std::to_string(cycle) + ".vtu");
data_out.write_vtu(output);
}
if (true) {
system_matrix.print(std::cout);
}
}
template <int dim> void Step6<dim>::run() {
for (unsigned int cycle = 0; cycle < 1; ++cycle) {
std::cout << "Cycle " << cycle << ':' << std::endl;
if (cycle == 0) {
Point<dim> p0;
p0(0) = -1;
p0(1) = -1;
// p0(2) = -1;
Point<dim> p1;
p1(0) = 1;
p1(1) = 1;
// p0(2) = 1;
const unsigned int divisions_integer = 2;
std::vector<unsigned int> divisions(dim);
divisions[0] = divisions_integer;
divisions[1] = divisions_integer;
GridGenerator::subdivided_hyper_rectangle(triangulation, divisions, p0,
p1);
triangulation.refine_global(1);
if (move_vertex) {
for (const auto &cell : triangulation.active_cell_iterators()) {
for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell;
++i) {
Point<dim> &v = cell->vertex(i);
if ((std::abs(v(0)) + std::abs(v(1))) < 1e-5)
v(1) -= displacement_vertex_a;
if ((std::abs(v(0) - 0.5) < 1e-5) && (std::abs(v(1)) < 1e-5))
v(1) -= displacement_vertex_b;
}
}
}
} else
refine_grid();
std::cout << " Number of active cells: "
<< triangulation.n_active_cells() << std::endl;
setup_system();
std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
assemble_system();
solve();
output_results(cycle);
}
}
int main() {
try {
Step6<2> laplace_problem_2d;
laplace_problem_2d.run();
} catch (std::exception &exc) {
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
} catch (...) {
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
