PolyFEM
Loading...
Searching...
No Matches
MVPolygonalBasis2d.cpp
Go to the documentation of this file.
1
2#include "MVPolygonalBasis2d.hpp"
3#include "BarycentricBasis2d.hpp"
4
5#include <memory>
6
7namespace polyfem
8{
9 namespace basis
10 {
11
12 namespace
13 {
14 template <typename Expr>
15 inline Eigen::RowVector2d rotatePi_2(const Expr &p) // rotation of pi/2
16 {
17 return Eigen::RowVector2d(-p(1), p(0));
18 }
19 } // anonymous namespace
20
22
23 void MVPolygonalBasis2d::meanvalue(const Eigen::MatrixXd &polygon, const Eigen::RowVector2d &point, Eigen::MatrixXd &b, const double tol)
24 {
25 const int n_boundary = polygon.rows();
26
27 Eigen::MatrixXd segments(n_boundary, 2);
28 Eigen::VectorXd radii(n_boundary);
29 Eigen::VectorXd areas(n_boundary);
30 Eigen::VectorXd products(n_boundary);
31 Eigen::VectorXd tangents(n_boundary);
32 Eigen::Matrix2d mat;
33
34 b.resize(n_boundary, 1);
35 b.setZero();
36
37 for (int i = 0; i < n_boundary; ++i)
38 {
39 segments.row(i) = polygon.row(i) - point;
40
41 radii(i) = segments.row(i).norm();
42
43 // we are on the vertex
44 if (radii(i) < tol)
45 {
46 b(i) = 1;
47
48 return;
49 }
50 }
51
52 for (int i = 0; i < n_boundary; ++i)
53 {
54 const int ip1 = (i + 1) == n_boundary ? 0 : (i + 1);
55
56 mat.row(0) = segments.row(i);
57 mat.row(1) = segments.row(ip1);
58
59 areas(i) = mat.determinant();
60 products(i) = segments.row(i).dot(segments.row(ip1));
61
62 // we are on the edge
63 if (fabs(areas[i]) < tol && products(i) < 0)
64 {
65 const double denominator = 1.0 / (radii(i) + radii(ip1));
66
67 b(i) = radii(ip1) * denominator;
68 b(ip1) = radii(i) * denominator;
69
70 return;
71 }
72 }
73
74 for (int i = 0; i < n_boundary; ++i)
75 {
76 const int ip1 = (i + 1) == n_boundary ? 0 : (i + 1);
77
78 tangents(i) = areas(i) / (radii(i) * radii(ip1) + products(i));
79 }
80
81 double W = 0;
82 for (int i = 0; i < n_boundary; ++i)
83 {
84 const int im1 = i == 0 ? (n_boundary - 1) : (i - 1);
85
86 b(i) = (tangents(im1) + tangents(i)) / radii(i);
87 W += b(i);
88 }
89
90 b /= W;
91 }
92
93 void MVPolygonalBasis2d::meanvalue_derivative(const Eigen::MatrixXd &polygon, const Eigen::RowVector2d &point, Eigen::MatrixXd &derivatives, const double tol)
94 {
95 const int n_boundary = polygon.rows();
96
97 // b.resize(n_boundary*n_points);
98 // std::fill(b.begin(), b.end(), 0);
99
100 derivatives.resize(n_boundary, 2);
101 derivatives.setZero();
102
103 Eigen::MatrixXd segments(n_boundary, 2);
104 Eigen::VectorXd radii(n_boundary);
105 Eigen::VectorXd areas(n_boundary);
106 Eigen::VectorXd products(n_boundary);
107 Eigen::VectorXd tangents(n_boundary);
108 Eigen::Matrix2d mat;
109
110 Eigen::MatrixXd areas_prime(n_boundary, 2);
111 Eigen::MatrixXd products_prime(n_boundary, 2);
112 Eigen::MatrixXd radii_prime(n_boundary, 2);
113 Eigen::MatrixXd tangents_prime(n_boundary, 2);
114 Eigen::MatrixXd w_prime(n_boundary, 2);
115
116 // Eigen::MatrixXd b(n_boundary, 1);
117
118 for (int i = 0; i < n_boundary; ++i)
119 {
120 segments.row(i) = polygon.row(i) - point;
121
122 radii(i) = segments.row(i).norm();
123
124 // we are on the vertex
125 if (radii(i) < tol)
126 {
127 // assert(false);
128 // b(i) = 1;
129 return;
130 }
131 }
132
133 int on_edge = -1;
134 double w0 = 0, w1 = 0;
135
136 for (int i = 0; i < n_boundary; ++i)
137 {
138 const int ip1 = (i + 1) == n_boundary ? 0 : (i + 1);
139
140 mat.row(0) = segments.row(i);
141 mat.row(1) = segments.row(ip1);
142
143 areas(i) = mat.determinant();
144 products(i) = segments.row(i).dot(segments.row(ip1));
145
146 // we are on the edge
147 if (fabs(areas[i]) < tol && products(i) < 0)
148 {
149 // const double denominator = 1.0/(radii(i) + radii(ip1));
150 // w0 = radii(ip1); // * denominator;
151 // w1 = radii(i); // * denominator;
152
153 // //https://link.springer.com/article/10.1007/s00371-013-0889-y
154 // const Eigen::RowVector2d NE = rotatePi_2(polygon.row(ip1) - polygon.row(i));
155 // const double sqrlengthE = NE.squaredNorm();
156
157 // const Eigen::RowVector2d N0 = rotatePi_2(point - polygon.row(i));
158 // const Eigen::RowVector2d N1 = rotatePi_2( polygon.row(ip1) - point);
159 // const double N0norm = N0.norm();
160 // const double N1norm = N1.norm();
161
162 // const Eigen::RowVector2d gradw0 = (N0.dot(N1) / (2*N0norm*N0norm*N0norm) + 1./(2.*N1norm) + 1./N0norm - 1./N1norm ) * NE / sqrlengthE;
163 // const Eigen::RowVector2d gradw1 = (1./(2*N1norm) + N0.dot(N1)/(2*N1norm*N1norm*N1norm) - 1./N0norm + 1./N1norm ) * NE / sqrlengthE;
164
165 // w_prime.setZero();
166 // w_prime.row(i) = gradw0;
167 // w_prime.row(ip1) = gradw1;
168
169 // assert(on_edge == -1);
170 // on_edge = i;
171 // continue;
172
173 // w_gradients_on_edges[e] = std::pair<point_t,point_t>(gradw0,gradw1);
174
175 // w_gradients[e0] += gradw0;
176 // w_gradients[e1] += gradw1;
177 // assert(false);
178 return;
179 }
180 const Eigen::RowVector2d vi = polygon.row(i);
181 const Eigen::RowVector2d vip1 = polygon.row(ip1);
182
183 areas_prime(i, 0) = vi(1) - vip1(1);
184 areas_prime(i, 1) = vip1(0) - vi(0);
185
186 products_prime.row(i) = 2 * point - vi - vip1;
187 radii_prime.row(i) = (point - vi) / radii(i);
188 }
189
190 for (int i = 0; i < n_boundary; ++i)
191 {
192 // if(i == on_edge)
193 // continue;
194
195 const int ip1 = (i + 1) == n_boundary ? 0 : (i + 1);
196
197 const double denominator = radii(i) * radii(ip1) + products(i);
198 const Eigen::RowVector2d denominator_prime = radii_prime.row(i) * radii(ip1) + radii(i) * radii_prime.row(ip1) + products_prime.row(i);
199
200 tangents_prime.row(i) = (areas_prime.row(i) * denominator - areas(i) * denominator_prime) / (denominator * denominator);
201 tangents(i) = areas(i) / denominator;
202 }
203
204 double W = 0;
205 Eigen::RowVector2d W_prime;
206 W_prime.setZero();
207
208 for (int i = 0; i < n_boundary; ++i)
209 {
210 const int im1 = (i > 0) ? (i - 1) : (n_boundary - 1);
211
212 if (i != on_edge && im1 != on_edge)
213 w_prime.row(i) = ((tangents_prime.row(im1) + tangents_prime.row(i)) * radii(i) - (tangents(im1) + tangents(i)) * radii_prime.row(i)) / (radii(i) * radii(i));
214 ;
215
216 W_prime += w_prime.row(i);
217 if (i == on_edge)
218 W += w0;
219 else if (im1 == on_edge)
220 W += w1;
221 else if (on_edge < 0)
222 W += (tangents(im1) + tangents(i)) / radii(i);
223 }
224
225 for (int i = 0; i < n_boundary; ++i)
226 {
227 const int im1 = (i > 0) ? (i - 1) : (n_boundary - 1);
228
229 double wi;
230 if (i == on_edge)
231 wi = w0;
232 else if (im1 == on_edge)
233 wi = w1;
234 else if (on_edge < 0)
235 wi = (tangents(im1) + tangents(i)) / radii(i);
236
237 derivatives.row(i) = (w_prime.row(i) * W - wi * W_prime) / (W * W);
238 }
239 }
240
241 int MVPolygonalBasis2d::build_bases(
242 const std::string &assembler_name,
243 const int dim,
244 const mesh::Mesh2D &mesh,
245 const int n_bases,
246 const int quadrature_order,
247 const int mass_quadrature_order,
248 std::vector<ElementBases> &bases,
249 std::vector<mesh::LocalBoundary> &local_boundary,
250 std::map<int, Eigen::MatrixXd> &mapped_boundary)
251 {
252 return BarycentricBasis2d::build_bases(assembler_name, dim, mesh, n_bases, quadrature_order, mass_quadrature_order, meanvalue, meanvalue_derivative, bases, local_boundary, mapped_boundary);
253 }
254
255 } // namespace basis
256} // namespace polyfem
mat
Eigen::MatrixXd mat
Definition TransientForm.cpp:24
polyfem::basis::BarycentricBasis2d::build_bases
static int build_bases(const std::string &assembler_name, const int dim, const mesh::Mesh2D &mesh, const int n_bases, const int quadrature_order, const int mass_quadrature_order, const std::function< void(const Eigen::MatrixXd &, const Eigen::RowVector2d &, Eigen::MatrixXd &, const double)> bc, const std::function< void(const Eigen::MatrixXd &, const Eigen::RowVector2d &, Eigen::MatrixXd &, const double)> bc_prime, std::vector< ElementBases > &bases, std::vector< mesh::LocalBoundary > &local_boundary, std::map< int, Eigen::MatrixXd > &mapped_boundary)
Definition BarycentricBasis2d.cpp:116
polyfem::basis::MVPolygonalBasis2d::build_bases
static int build_bases(const std::string &assembler_name, const int dim, const mesh::Mesh2D &mesh, const int n_bases, const int quadrature_order, const int mass_quadrature_order, std::vector< ElementBases > &bases, std::vector< mesh::LocalBoundary > &local_boundary, std::map< int, Eigen::MatrixXd > &mapped_boundary)
Definition MVPolygonalBasis2d.cpp:241
polyfem::basis::MVPolygonalBasis2d::meanvalue
static void meanvalue(const Eigen::MatrixXd &polygon, const Eigen::RowVector2d &point, Eigen::MatrixXd &b, const double tol)
Definition MVPolygonalBasis2d.cpp:23
polyfem::basis::MVPolygonalBasis2d::meanvalue_derivative
static void meanvalue_derivative(const Eigen::MatrixXd &polygon, const Eigen::RowVector2d &point, Eigen::MatrixXd &derivatives, const double tol)
Definition MVPolygonalBasis2d.cpp:93