34 #include "DGtal/base/Common.h"
35 #include "ConfigTest.h"
36 #include "DGtalCatch.h"
37 #include "DGtal/helpers/StdDefs.h"
39 #include "DGtal/dec/PolygonalCalculus.h"
40 #include "DGtal/shapes/SurfaceMesh.h"
41 #include "DGtal/shapes/MeshHelpers.h"
42 #include "DGtal/helpers/Shortcuts.h"
43 #include "DGtal/math/linalg/DirichletConditions.h"
47 using namespace DGtal;
62 std::vector<RealPoint> positions = {
RealPoint( 0, 0, 0 ) ,
72 std::vector<Mesh::Vertices> faces = { { 1, 0, 2, 3 },
79 Mesh box(positions.cbegin(), positions.cend(),
80 faces.cbegin(), faces.cend());
84 SECTION(
"Construction and basic operators")
90 auto x = boxCalculus.
X(f);
91 auto d = boxCalculus.
D(f);
92 auto a = boxCalculus.
A(f);
118 for(
auto f=0; f<6; ++f)
121 box.computeFaceNormalsFromPositions();
122 for(
auto f=0; f<6; ++f)
125 auto n = box.faceNormal(f);
137 auto d = boxCalculus.
D(f);
141 phi << 1.0, 3.0, 2.0, 6.0;
142 expected << 2,-1,4,-5;
148 SECTION(
"Structural propertes")
153 phi << 1.0, 3.0, 2.0, 6.0;
158 auto cogphi = coG*phi;
161 REQUIRE( gphi.dot(cogphi) == 0.0);
171 auto P = boxCalculus.
P(f);
181 auto curl = boxCalculus.
curl(f);
186 SECTION(
"Local Laplace-Beltrami")
189 auto nf = box.incidentVertices(f).size();
193 phi << 1.0, 1.0, 1.0, 1.0;
198 SECTION(
"Check lumped mass matrix")
202 for(
auto v=0; v < box.nbVertices(); ++v )
203 a += M.coeffRef(v,v);
206 for(
auto f=0; f < box.
nbFaces(); ++f )
207 fa += box.faceArea(f);
223 trace.
info()<< boxCalculusCached <<std::endl;
227 for(
auto i=0; i < 1000 ; ++i)
233 for(
auto i=0; i < 1000 ; ++i)
237 REQUIRE(L.norm() == Approx(LC.norm()));
243 TEST_CASE(
"Testing PolygonalCalculus and DirichletConditions" )
254 params(
"polynomial",
"0.1*y*y -0.1*x*x - 2.0*z" )(
"gridstep", 2.0 );
262 std::vector<std::vector< Index > > faces;
263 std::vector<RealPoint> positions = primalSurface->positions();
264 for(
auto face= 0 ; face < primalSurface->nbFaces(); ++face)
265 faces.push_back(primalSurface->incidentVertices( face ));
267 Mesh surfmesh =
Mesh( positions.begin(), positions.end(),
268 faces.begin(), faces.end() );
275 REQUIRE( boundaryEdges.size() == 168 );
279 PolyDEC calculus( surfmesh );
285 DC::IntegerVector b = DC::IntegerVector::Zero( g.rows() );
287 SECTION(
"Solve Poisson problem with boundary Dirichlet conditions")
289 for (
double scale = 0.1; scale < 2.0; scale *= 2.0 )
291 std::cout <<
"scale=" << scale << std::endl;
292 auto phi = [&](
Index v)
294 return cos(scale*(surfmesh.
position(v)[0]))
298 for(
auto &e: boundaryEdges)
301 auto v1 = adjVertices.first;
302 auto v2 = adjVertices.second;
309 PolyDEC::Solver solver;
311 solver.compute( L_dirichlet );
312 REQUIRE( solver.info() == Eigen::Success );
315 REQUIRE( solver.info() == Eigen::Success );
317 double min_phi = 0.0;
318 double max_phi = 0.0;
321 double min_i_u = 0.0;
322 double max_i_u = 0.0;
325 for (
auto v = 0; v < surfmesh.
nbVertices(); ++v )
327 min_phi = std::min( min_phi, phi( v ) );
328 max_phi =
std::max( max_phi, phi( v ) );
329 min_u = std::min( min_u , u ( v ) );
330 max_u =
std::max( max_u , u ( v ) );
333 min_i_u = std::min( min_i_u, u ( v ) );
334 max_i_u =
std::max( max_i_u, u ( v ) );