Skip to content

Commit 8bd1bfd

Browse files
author
Teseo Schneider
committed
qp done
1 parent 76ccb1b commit 8bd1bfd

2 files changed

Lines changed: 154 additions & 64 deletions

File tree

‎src/min_quad_with_fixed.cpp‎

Lines changed: 83 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,83 @@
1+
//TODO: __example remove __copy
2+
3+
#include <common.h>
4+
#include <npe.h>
5+
#include <typedefs.h>
6+
7+
8+
9+
10+
11+
12+
#include <igl/min_quad_with_fixed.h>
13+
14+
const char *ds_min_quad_with_fixed = R"igl_Qu8mg5v7(
15+
16+
MIN_QUAD_WITH_FIXED Minimize a quadratic energy of the form
17+
trace( 0.5*Z'*A*Z + Z'*B + constant )
18+
subject to
19+
Z(known,:) = Y, and
20+
Aeq*Z = Beq
21+
22+
Parameters
23+
----------
24+
A n by n matrix of quadratic coefficients
25+
B n by 1 column of linear coefficients
26+
known list of indices to known rows in Z
27+
Y list of fixed values corresponding to known rows in Z
28+
Aeq m by n list of linear equality constraint coefficients
29+
Beq m by 1 list of linear equality constraint constant values
30+
is_A_pd flag specifying whether A(unknown,unknown) is positive definite
31+
32+
33+
Returns
34+
-------
35+
Z n by k solution
36+
37+
See also
38+
--------
39+
40+
41+
Notes
42+
-----
43+
None
44+
45+
Examples
46+
--------
47+
48+
49+
)igl_Qu8mg5v7";
50+
51+
npe_function(min_quad_with_fixed)
52+
npe_doc(ds_min_quad_with_fixed)
53+
54+
npe_arg(A, sparse_double)
55+
npe_arg(B, dense_double)
56+
npe_arg(known, dense_int, dense_long, dense_longlong)
57+
npe_arg(Y, dense_double)
58+
npe_arg(Aeq, sparse_double)
59+
npe_arg(Beq, dense_double)
60+
npe_arg(is_A_pd, bool)
61+
62+
63+
npe_begin_code()
64+
65+
assert_nonzero_rows(A, "A");
66+
if(Aeq.size() > 0)
67+
assert_cols_match(A, Aeq, "A", "Aeq");
68+
assert_rows_match(A, B, "A", "B");
69+
assert_cols_match(B, Y, "B", "Y");
70+
if (Beq.size() > 0)
71+
assert_cols_match(B, Beq, "B", "Beq");
72+
assert_rows_match(Aeq, Beq, "Aeq", "Beq");
73+
74+
Eigen::SparseMatrix<double> A_copy = A;
75+
Eigen::SparseMatrix<double> Aeq_copy = Aeq;
76+
77+
Eigen::MatrixXd sol;
78+
79+
bool ok = igl::min_quad_with_fixed(A_copy, B, known, Y, Aeq_copy, Beq, is_A_pd, sol);
80+
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> sol_row_major = sol;
81+
return std::make_pair(ok, npe::move(sol_row_major));
82+
83+
npe_end_code()

‎tutorial/tutorials.ipynb‎

Lines changed: 71 additions & 64 deletions
Original file line numberDiff line numberDiff line change
@@ -825,38 +825,49 @@
825825
"metadata": {},
826826
"outputs": [],
827827
"source": [
828-
"# v, f = igl.read_triangle_mesh(\"data/cheburashka.off\")\n",
829-
"\n",
830-
"# ## Two fixed points: Left hand, left foot should have values 1 and -1\n",
831-
"# b = np.array([4331, 5957])\n",
832-
"# bc = np.array([1, -1])\n",
833-
"\n",
834-
"# ## Construct Laplacian and mass matrix\n",
835-
"# l = igl.cotmatrix(v, f)\n",
836-
"# m = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_VORONOI)\n",
837-
"# minv = sp.sparse.diags(1 / m.diagonal())\n",
838-
"\n",
839-
"# ## Bi-Laplacian\n",
840-
"# q = l.dot(minv.dot(l))\n",
841-
"\n",
842-
"# ## Solve with only equality constraints\n",
843-
"# z1 = igl.min_quad_with_fixed(q, b, bc)\n",
844-
"\n",
845-
"# ## Solve with equality and linear constraints\n",
846-
"# bl = np.array([[6074, 6523]])\n",
847-
"# z2 = igl.min_quad_with_fixed(q, b, bc, bl)\n",
848-
"\n",
849-
"# ## Normalize colors to same range\n",
850-
"# min_z = min(np.min(z1), np.min(z2))\n",
851-
"# max_z = max(np.max(z1), np.max(z2))\n",
852-
"# z = [(z1 - min_z) / (max_z - min_z), (z2 - min_z) / (max_z - min_z)]\n",
853-
" \n",
854-
"# ## Plot the functions\n",
855-
"# p = plot(v, f, z1, shading={\"wireframe\":False}, return_plot=True)\n",
856-
"\n",
857-
"# @interact(function=[('z0', 0), ('z1', 1)])\n",
858-
"# def sf(function):\n",
859-
"# p.update_object(colors=z[function])"
828+
"import igl\n",
829+
"import numpy as np\n",
830+
"import scipy as sp\n",
831+
"\n",
832+
"v, f = igl.read_triangle_mesh(\"data/cheburashka.off\")\n",
833+
"\n",
834+
"## Two fixed points: Left hand, left foot should have values 1 and -1\n",
835+
"b = np.array([4331, 5957])\n",
836+
"bc = np.array([1., -1.])\n",
837+
"B = np.zeros((v.shape[0], 1))\n",
838+
"\n",
839+
"## Construct Laplacian and mass matrix\n",
840+
"L = igl.cotmatrix(v, f)\n",
841+
"M = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_VORONOI)\n",
842+
"Minv = sp.sparse.diags(1 / M.diagonal())\n",
843+
"\n",
844+
"## Bi-Laplacian\n",
845+
"Q = L @ (Minv @ L)\n",
846+
"Q.sort_indices()\n",
847+
"\n",
848+
"## Solve with only equality constraints\n",
849+
"Aeq = sp.sparse.csc_matrix((0, 0))\n",
850+
"Beq = np.array([])\n",
851+
"_, z1 = igl.min_quad_with_fixed(Q, B, b, bc, Aeq, Beq, True)\n",
852+
"\n",
853+
"## Solve with equality and linear constraints\n",
854+
"Aeq = sp.sparse.csc_matrix((1, v.shape[0]))\n",
855+
"Aeq[0,6074] = 1\n",
856+
"Aeq[0, 6523] = -1\n",
857+
"Beq = np.array([0.])\n",
858+
"_, z2 = igl.min_quad_with_fixed(Q, B, b, bc, Aeq, Beq, True)\n",
859+
"\n",
860+
"## Normalize colors to same range\n",
861+
"min_z = min(np.min(z1), np.min(z2))\n",
862+
"max_z = max(np.max(z1), np.max(z2))\n",
863+
"z = [(z1 - min_z) / (max_z - min_z), (z2 - min_z) / (max_z - min_z)]\n",
864+
"\n",
865+
"## Plot the functions\n",
866+
"p = plot(v, f, z1, shading={\"wireframe\":False}, return_plot=True)\n",
867+
"\n",
868+
"@interact(function=[('z0', 0), ('z1', 1)])\n",
869+
"def sf(function):\n",
870+
" p.update_object(colors=z[function])"
860871
]
861872
},
862873
{
@@ -1120,7 +1131,7 @@
11201131
"metadata": {},
11211132
"outputs": [],
11221133
"source": [
1123-
"v, _, _, f, _, _ = igl.read_obj(\"data/bump-domain.obj\")\n",
1134+
"v, f = igl.read_triangle_mesh(\"data/bump-domain.obj\")\n",
11241135
"u = v.copy()\n",
11251136
"\n",
11261137
"# Find boundary vertices outside annulus\n",
@@ -1262,7 +1273,7 @@
12621273
"metadata": {},
12631274
"outputs": [],
12641275
"source": [
1265-
"v, f, _ = igl.read_off(\"data/camelhead.off\")\n",
1276+
"v, f = igl.read_triangle_mesh(\"data/camelhead.off\")\n",
12661277
"\n",
12671278
"# Fix two points on the boundary\n",
12681279
"b = np.array([2, 1])\n",
@@ -1418,12 +1429,12 @@
14181429
"metadata": {},
14191430
"outputs": [],
14201431
"source": [
1421-
"v, f, _ = igl.read_off(\"data/fertility.off\")\n",
1432+
"v, f = igl.read_triangle_mesh(\"data/fertility.off\")\n",
14221433
"\n",
14231434
"n = igl.per_vertex_normals(v, f)\n",
14241435
"\n",
14251436
"# Compute ambient occlusion factor using embree\n",
1426-
"ao = igl.ambient_occlusion(v, f, v, n, 500)\n",
1437+
"ao = igl.ambient_occlusion(v, f, v, n, 50)\n",
14271438
"ao = 1.0 - ao\n",
14281439
"\n",
14291440
"plot(v, f, ao, shading={\"colormap\": \"gist_gray\"})"
@@ -1453,33 +1464,7 @@
14531464
"\n",
14541465
"The example computes these quantities and\n",
14551466
"does a basic statistic analysis that allows to estimate the isometry and\n",
1456-
"regularity of a mesh:\n",
1457-
"\n",
1458-
"```bash\n",
1459-
"Irregular vertices:\n",
1460-
"136/2400 (5.67%)\n",
1461-
"Areas (Min/Max)/Avg_Area Sigma:\n",
1462-
"0.01/5.33 (0.87)\n",
1463-
"Angles in degrees (Min/Max) Sigma:\n",
1464-
"17.21/171.79 (15.36)\n",
1465-
"```\n",
1466-
"\n",
1467-
"The first row contains the number and percentage of irregular vertices, which\n",
1468-
"is particularly important for quadrilateral meshes when they are used to define\n",
1469-
"subdivision surfaces: every singular point will result in a point of the\n",
1470-
"surface that is only C^1.\n",
1471-
"\n",
1472-
"The second row reports the area of the minimal element, maximal element and the\n",
1473-
"standard deviation. These numbers are normalized by the mean area, so in the\n",
1474-
"example above 5.33 max area means that the biggest face is 5 times larger than\n",
1475-
"the average face. An ideal isotropic mesh would have both min and max area\n",
1476-
"close to 1.\n",
1477-
"\n",
1478-
"The third row measures the face angles, which should be close to 60 degrees (90\n",
1479-
"for quads) in a perfectly regular triangulation. For FEM purposes, the closer\n",
1480-
"the angles are to 60 degrees the more stable will the optimization be. In this\n",
1481-
"case, it is clear that the mesh is of bad quality and it will probably result\n",
1482-
"in artifacts if used for solving PDEs."
1467+
"regularity of a mesh:"
14831468
]
14841469
},
14851470
{
@@ -1522,6 +1507,28 @@
15221507
"print(\"Angles in degrees (Min/Max) Sigma: \\n%.2f/%.2f (%.2f)\\n\"%(angle_min, angle_max, angle_sigma))"
15231508
]
15241509
},
1510+
{
1511+
"cell_type": "markdown",
1512+
"metadata": {},
1513+
"source": [
1514+
"The first row contains the number and percentage of irregular vertices, which\n",
1515+
"is particularly important for quadrilateral meshes when they are used to define\n",
1516+
"subdivision surfaces: every singular point will result in a point of the\n",
1517+
"surface that is only C^1.\n",
1518+
"\n",
1519+
"The second row reports the area of the minimal element, maximal element and the\n",
1520+
"standard deviation. These numbers are normalized by the mean area, so in the\n",
1521+
"example above 5.33 max area means that the biggest face is 5 times larger than\n",
1522+
"the average face. An ideal isotropic mesh would have both min and max area\n",
1523+
"close to 1.\n",
1524+
"\n",
1525+
"The third row measures the face angles, which should be close to 60 degrees (90\n",
1526+
"for quads) in a perfectly regular triangulation. For FEM purposes, the closer\n",
1527+
"the angles are to 60 degrees the more stable will the optimization be. In this\n",
1528+
"case, it is clear that the mesh is of bad quality and it will probably result\n",
1529+
"in artifacts if used for solving PDEs."
1530+
]
1531+
},
15251532
{
15261533
"cell_type": "markdown",
15271534
"metadata": {},
@@ -1617,7 +1624,7 @@
16171624
"name": "python",
16181625
"nbconvert_exporter": "python",
16191626
"pygments_lexer": "ipython3",
1620-
"version": "3.7.3"
1627+
"version": "3.6.7"
16211628
}
16221629
},
16231630
"nbformat": 4,

0 commit comments

Comments
 (0)