Wildmeshing Toolkit
Loading...
Searching...
No Matches
periodic_optimization.cpp
Go to the documentation of this file.
1#include "periodic_optimization.hpp"
2
3#include <wmtk/Mesh.hpp>
4#include <wmtk/Scheduler.hpp>
5#include <wmtk/TetMesh.hpp>
6#include <wmtk/TriMesh.hpp>
7
8#include <wmtk/components/utils/get_attributes.hpp>
9#include <wmtk/multimesh/consolidate.hpp>
10#include <wmtk/utils/Logger.hpp>
11
12
13#include <wmtk/operations/attribute_new/SplitNewAttributeStrategy.hpp>
14#include <wmtk/operations/attribute_update/AttributeTransferStrategy.hpp>
15
16#include <wmtk/operations/AMIPSOptimizationSmoothing.hpp>
17#include <wmtk/operations/AMIPSOptimizationSmoothingPeriodic.hpp>
18#include <wmtk/operations/AndOperationSequence.hpp>
19#include <wmtk/operations/EdgeCollapse.hpp>
20#include <wmtk/operations/EdgeSplit.hpp>
21#include <wmtk/operations/MinOperationSequence.hpp>
22#include <wmtk/operations/OptimizationSmoothing.hpp>
23#include <wmtk/operations/OrOperationSequence.hpp>
24#include <wmtk/operations/Rounding.hpp>
25#include <wmtk/operations/TetWildTangentialLaplacianSmoothing.hpp>
26#include <wmtk/operations/composite/ProjectOperation.hpp>
27#include <wmtk/operations/composite/TetEdgeSwap.hpp>
28#include <wmtk/operations/composite/TetFaceSwap.hpp>
29#include <wmtk/operations/composite/TriEdgeSwap.hpp>
30
31
32#include <wmtk/function/LocalNeighborsSumFunction.hpp>
33#include <wmtk/function/PerSimplexFunction.hpp>
34#include <wmtk/function/simplex/AMIPS.hpp>
35
36#include <wmtk/function/utils/amips.hpp>
37#include <wmtk/invariants/CollapseEnergyBeforeInvariantDouble.hpp>
38#include <wmtk/invariants/EdgeValenceInvariant.hpp>
39#include <wmtk/invariants/EnergyFilterInvariant.hpp>
40#include <wmtk/invariants/EnvelopeInvariant.hpp>
41#include <wmtk/invariants/FunctionInvariant.hpp>
42#include <wmtk/invariants/InteriorEdgeInvariant.hpp>
43#include <wmtk/invariants/InteriorSimplexInvariant.hpp>
44#include <wmtk/invariants/InteriorVertexInvariant.hpp>
45#include <wmtk/invariants/MaxFunctionInvariant.hpp>
46#include <wmtk/invariants/MultiMeshLinkConditionInvariant.hpp>
47#include <wmtk/invariants/NoBoundaryCollapseToInteriorInvariant.hpp>
48#include <wmtk/invariants/NoChildMeshAttachingInvariant.hpp>
49#include <wmtk/invariants/RoundedInvariant.hpp>
50#include <wmtk/invariants/SeparateSubstructuresInvariant.hpp>
51#include <wmtk/invariants/SimplexInversionInvariant.hpp>
52#include <wmtk/invariants/Swap32EnergyBeforeInvariantDouble.hpp>
53#include <wmtk/invariants/Swap44EnergyBeforeInvariantDouble.hpp>
54#include <wmtk/invariants/Swap56EnergyBeforeInvariantDouble.hpp>
55#include <wmtk/invariants/TodoInvariant.hpp>
56
57#include <wmtk/multimesh/utils/extract_child_mesh_from_tag.hpp>
58
59#include <wmtk/utils/orient.hpp>
60
61#include <wmtk/io/MeshReader.hpp>
62#include <wmtk/io/ParaviewWriter.hpp>
63
64#include <queue>
65#include <wmtk/simplex/k_ring.hpp>
66#include <wmtk/simplex/link.hpp>
67
68#include <fstream>
69
70
71namespace wmtk::components::internal {
72
73using namespace simplex;
74using namespace operations;
75using namespace operations::tri_mesh;
76using namespace operations::tet_mesh;
77using namespace operations::composite;
78using namespace function;
79using namespace invariants;
80
81void write(
82 const Mesh& mesh,
83 const std::string& out_dir,
84 const std::string& name,
85 const std::string& vname,
86 const int64_t index,
87 const bool intermediate_output);
88
89void adjust_sizing_field(
90 Mesh& m,
91 const TypedAttributeHandle<double>& coordinate_handle,
92 const TypedAttributeHandle<double>& edge_length_handle,
93 const TypedAttributeHandle<double>& sizing_field_scalar_handle,
94 const TypedAttributeHandle<double>& energy_handle,
95 const TypedAttributeHandle<double>& target_edge_length_handle,
96 const TypedAttributeHandle<char>& visited_handle,
97 const double stop_energy,
98 const double current_max_energy,
99 const double initial_target_edge_length,
100 const double min_target_edge_length);
101
102void set_operation_energy_filter(
103 Mesh& m,
104 const TypedAttributeHandle<double>& coordinate_handle,
105 const TypedAttributeHandle<double>& energy_handle,
106 const TypedAttributeHandle<char>& energy_filter_handle,
107 const TypedAttributeHandle<char>& visited_handle,
108 const double stop_energy,
109 const double current_max_energy,
110 const double initial_target_edge_length);
111
112void periodic_optimization(
113 Mesh& periodic_mesh,
114 Mesh& position_mesh,
115 Mesh& surface_mesh,
116 const double target_edge_length,
117 const double max_amips_energy,
118 const int64_t passes,
119 const double envelope_size,
120 const bool intermediate_output,
121 std::vector<attribute::MeshAttributeHandle>& pass_through_attributes,
122 std::string output_dir,
123 std::string output)
124{
125 auto position_handle =
126 position_mesh.get_attribute_handle<double>("vertices", PrimitiveType::Vertex);
127 auto surface_position_handle =
128 surface_mesh.get_attribute_handle<double>("vertices", PrimitiveType::Vertex);
129
130 auto position_accessor = position_mesh.create_const_accessor(position_handle.as<double>());
131
132 wmtk::logger().info(
133 "periodic mesh #v: {}",
134 periodic_mesh.get_all(PrimitiveType::Vertex).size());
135 wmtk::logger().info(
136 "position mesh #v: {}",
137 position_mesh.get_all(PrimitiveType::Vertex).size());
138
139 Eigen::Vector3d bmin, bmax;
140 bmin.setConstant(std::numeric_limits<double>::max());
141 bmax.setConstant(std::numeric_limits<double>::lowest());
142
143 const auto vertices = position_mesh.get_all(PrimitiveType::Vertex);
144 for (const auto& v : vertices) {
145 const auto p = position_accessor.vector_attribute(v);
146 for (int64_t d = 0; d < bmax.size(); ++d) {
147 bmin[d] = std::min(bmin[d], p[d]);
148 bmax[d] = std::max(bmax[d], p[d]);
149 }
150 }
151
152 const double bbdiag = (bmax - bmin).norm();
153 const double target_edge_length_abs = bbdiag * target_edge_length;
154
155 const double period_x = bmax[0] - bmin[0];
156 const double period_y = bmax[1] - bmin[1];
157 const double period_z = bmax[2] - bmin[2];
158
159
161 // envelope
163
164 auto envelope_invariant = std::make_shared<EnvelopeInvariant>(
165 surface_position_handle,
166 std::sqrt(2) * envelope_size,
167 surface_position_handle);
168
169 auto propagate_to_child_position =
170 [](const Eigen::MatrixX<double>& P) -> Eigen::VectorX<double> { return P; };
171
172 auto propagate_to_parent_position =
173 [](const Eigen::MatrixX<double>& P) -> Eigen::VectorX<double> {
174 assert(P.cols() == 1);
175 return P.col(0);
176 };
177
178 auto update_parent_position = std::make_shared<SingleAttributeTransferStrategy<double, double>>(
179 position_handle,
180 surface_position_handle,
181 propagate_to_parent_position);
182
183 auto update_child_position = std::make_shared<SingleAttributeTransferStrategy<double, double>>(
184 surface_position_handle,
185 position_handle,
186 propagate_to_child_position);
187
189 // amips energy
191
192 auto amips_handle =
193 position_mesh.register_attribute<double>("amips", PrimitiveType::Tetrahedron, 1);
194 auto amips_accessor = position_mesh.create_accessor(amips_handle.as<double>());
195
196 auto compute_amips = [](const Eigen::MatrixXd& P) -> Eigen::VectorXd {
197 // tet
198 std::array<double, 12> pts;
199 for (size_t i = 0; i < 4; ++i) {
200 for (size_t j = 0; j < 3; ++j) {
201 pts[3 * i + j] = P(j, i);
202 }
203 }
204 const double a = Tet_AMIPS_energy(pts);
205 return Eigen::VectorXd::Constant(1, a);
206 };
207 auto amips_update =
208 std::make_shared<wmtk::operations::SingleAttributeTransferStrategy<double, double>>(
209 amips_handle,
210 position_handle,
211 compute_amips);
212 amips_update->run_on_all();
213
214 double max_amips = std::numeric_limits<double>::lowest();
215 double min_amips = std::numeric_limits<double>::max();
216
217 for (const auto& t : position_mesh.get_all(PrimitiveType::Tetrahedron)) {
218 double e = amips_accessor.scalar_attribute(t);
219 max_amips = std::max(max_amips, e);
220 min_amips = std::min(min_amips, e);
221 }
222
223 logger().info("Initial Max AMIPS Energy: {}, Min AMIPS Energy: {}", max_amips, min_amips);
224
226 // target edge length
228
229 auto target_edge_length_handle = position_mesh.register_attribute<double>(
230 "target_edge_length",
231 PrimitiveType::Edge,
232 1,
233 false,
234 target_edge_length_abs); // defaults to target edge length
235
237 // edge length
239
240 auto edge_length_handle =
241 position_mesh.register_attribute<double>("edge_length", PrimitiveType::Edge, 1);
242 auto edge_length_accessor = position_mesh.create_accessor(edge_length_handle.as<double>());
243 // Edge length update
244 auto compute_edge_length = [](const Eigen::MatrixXd& P) -> Eigen::VectorXd {
245 return Eigen::VectorXd::Constant(1, sqrt((P.col(0) - P.col(1)).squaredNorm()));
246 };
247 auto edge_length_update =
248 std::make_shared<wmtk::operations::SingleAttributeTransferStrategy<double, double>>(
249 edge_length_handle,
250 position_handle,
251 compute_edge_length);
252 edge_length_update->run_on_all();
253
255 // sizing field scalar
257
258 auto sizing_field_scalar_handle = position_mesh.register_attribute<double>(
259 "sizing_field_scalar",
260 PrimitiveType::Vertex,
261 1,
262 false,
263 1); // defaults to 1
264
266 // sizing field update flags
268 auto visited_vertex_flag_handle =
269 position_mesh
270 .register_attribute<char>("visited_vertex", PrimitiveType::Vertex, 1, false, char(1));
271
273 // energy filter flag
275 auto energy_filter_handle =
276 position_mesh
277 .register_attribute<char>("energy_filter", PrimitiveType::Vertex, 1, false, char(1));
278
279 auto energy_filter_accessor = position_mesh.create_accessor<char>(energy_filter_handle);
280
281 auto update_energy_filter_func = [](const Eigen::MatrixX<double>& P) -> Eigen::VectorX<char> {
282 return Eigen::VectorX<char>::Constant(1, char(1));
283 };
284 auto energy_filter_update =
285 std::make_shared<wmtk::operations::SingleAttributeTransferStrategy<char, double>>(
286 energy_filter_handle,
287 position_handle,
288 update_energy_filter_func);
289
291 // renew flags
293 auto visited_edge_flag_handle =
294 position_mesh
295 .register_attribute<char>("visited_edge", PrimitiveType::Edge, 1, false, char(1));
296
297 auto update_flag_func = [](const Eigen::MatrixX<double>& P) -> Eigen::VectorX<char> {
298 return Eigen::VectorX<char>::Constant(1, char(1));
299 };
300 auto tag_update =
301 std::make_shared<wmtk::operations::SingleAttributeTransferStrategy<char, double>>(
302 visited_edge_flag_handle,
303 position_handle,
304 update_flag_func);
305
307 // pass through
309 pass_through_attributes.push_back(edge_length_handle);
310 pass_through_attributes.push_back(amips_handle);
311 pass_through_attributes.push_back(visited_vertex_flag_handle);
312 pass_through_attributes.push_back(energy_filter_handle);
313 pass_through_attributes.push_back(surface_position_handle);
314
316 // invariants
318
319 auto inversion_invariant = std::make_shared<SimplexInversionInvariant<double>>(
320 position_mesh,
321 position_handle.as<double>());
322
323 std::shared_ptr<function::PerSimplexFunction> amips =
324 std::make_shared<AMIPS>(position_mesh, position_handle);
325
326 auto function_invariant =
327 std::make_shared<MaxFunctionInvariant>(position_mesh.top_simplex_type(), amips);
328
329 auto link_condition = std::make_shared<MultiMeshLinkConditionInvariant>(position_mesh);
330
331 auto todo_larger = std::make_shared<TodoLargerInvariant>(
332 position_mesh,
333 edge_length_handle.as<double>(),
334 target_edge_length_handle.as<double>(),
335 4.0 / 3.0);
336
337 auto todo_smaller = std::make_shared<TodoSmallerInvariant>(
338 position_mesh,
339 edge_length_handle.as<double>(),
340 target_edge_length_handle.as<double>(),
341 4.0 / 5.0);
342
343 auto valence_3 = std::make_shared<EdgeValenceInvariant>(position_mesh, 3);
344 auto valence_4 = std::make_shared<EdgeValenceInvariant>(position_mesh, 4);
345 auto valence_5 = std::make_shared<EdgeValenceInvariant>(position_mesh, 5);
346
347 auto invariant_separate_substructures =
348 std::make_shared<invariants::SeparateSubstructuresInvariant>(position_mesh);
349
351 // operations
353
354 std::vector<std::shared_ptr<Operation>> ops;
355 std::vector<std::string> ops_name;
356
357 auto long_edges_first = [&](const simplex::Simplex& s) {
358 assert(s.primitive_type() == PrimitiveType::Edge);
359 return -edge_length_accessor.scalar_attribute(s.tuple());
360 };
361 auto short_edges_first = [&](const simplex::Simplex& s) {
362 assert(s.primitive_type() == PrimitiveType::Edge);
363 return edge_length_accessor.scalar_attribute(s.tuple());
364 };
365
367 // 1) EdgeSplit
369 auto split = std::make_shared<EdgeSplit>(position_mesh);
370 split->set_priority(long_edges_first);
371
372 split->add_invariant(todo_larger);
373 split->add_invariant(
374 std::make_shared<EnergyFilterInvariant>(position_mesh, energy_filter_handle.as<char>()));
375 split->add_invariant(inversion_invariant);
376
377 split->set_new_attribute_strategy(position_handle);
378 split->set_new_attribute_strategy(sizing_field_scalar_handle);
379 split->set_new_attribute_strategy(
380 visited_edge_flag_handle,
381 wmtk::operations::SplitBasicStrategy::None,
382 wmtk::operations::SplitRibBasicStrategy::None);
383 split->set_new_attribute_strategy(
384 target_edge_length_handle,
385 wmtk::operations::SplitBasicStrategy::Copy,
386 wmtk::operations::SplitRibBasicStrategy::Mean);
387
388 for (const auto& attr : pass_through_attributes) {
389 split->set_new_attribute_strategy(
390 attr,
391 wmtk::operations::SplitBasicStrategy::None,
392 wmtk::operations::SplitRibBasicStrategy::None);
393 }
394
395 split->add_transfer_strategy(amips_update);
396 split->add_transfer_strategy(edge_length_update);
397 split->add_transfer_strategy(tag_update); // for renew the queue
398 split->add_transfer_strategy(energy_filter_update);
399 split->add_transfer_strategy(update_child_position);
400
401 ops.emplace_back(split);
402 ops_name.emplace_back("SPLIT");
403
405 // 1) EdgeCollapse
407
408 auto clps_strat1 = std::make_shared<CollapseNewAttributeStrategy<double>>(position_handle);
409 clps_strat1->set_strategy(CollapseBasicStrategy::CopyOther);
410
411 auto clps_strat2 = std::make_shared<CollapseNewAttributeStrategy<double>>(position_handle);
412 clps_strat2->set_strategy(CollapseBasicStrategy::CopyTuple);
413
414 auto setup_collapse = [&](std::shared_ptr<EdgeCollapse>& collapse) {
415 collapse->add_invariant(invariant_separate_substructures);
416 collapse->add_invariant(link_condition);
417 collapse->add_invariant(inversion_invariant);
418 collapse->add_invariant(envelope_invariant);
419
420 collapse->set_new_attribute_strategy(
421 visited_edge_flag_handle,
422 wmtk::operations::CollapseBasicStrategy::None);
423
424 collapse->add_transfer_strategy(tag_update);
425 collapse->add_transfer_strategy(energy_filter_update);
426 for (const auto& attr : pass_through_attributes) {
427 collapse->set_new_attribute_strategy(
428 attr,
429 wmtk::operations::CollapseBasicStrategy::None);
430 }
431 collapse->set_new_attribute_strategy(
432 target_edge_length_handle,
433 wmtk::operations::CollapseBasicStrategy::None);
434
435 collapse->add_transfer_strategy(amips_update);
436 collapse->add_transfer_strategy(edge_length_update);
437
438 collapse->add_transfer_strategy(update_child_position);
439 };
440
441 auto collapse1 = std::make_shared<EdgeCollapse>(position_mesh);
442 collapse1->add_invariant(std::make_shared<CollapseEnergyBeforeInvariantDouble>(
443 position_mesh,
444 position_handle.as<double>(),
445 amips_handle.as<double>(),
446 1));
447
448 collapse1->set_new_attribute_strategy(position_handle, clps_strat1);
449 collapse1->set_new_attribute_strategy(sizing_field_scalar_handle, clps_strat1);
450 setup_collapse(collapse1);
451
452 auto collapse2 = std::make_shared<EdgeCollapse>(position_mesh);
453 collapse2->add_invariant(std::make_shared<CollapseEnergyBeforeInvariantDouble>(
454 position_mesh,
455 position_handle.as<double>(),
456 amips_handle.as<double>(),
457 0));
458
459 collapse2->set_new_attribute_strategy(position_handle, clps_strat2);
460 collapse2->set_new_attribute_strategy(sizing_field_scalar_handle, clps_strat2);
461 setup_collapse(collapse2);
462
463 auto collapse = std::make_shared<OrOperationSequence>(position_mesh);
464 collapse->add_operation(collapse1);
465 collapse->add_operation(collapse2);
466 collapse->set_priority(short_edges_first);
467 collapse->add_invariant(
468 std::make_shared<EnergyFilterInvariant>(position_mesh, energy_filter_handle.as<char>()));
469 collapse->add_invariant(todo_smaller);
470
471 collapse->add_transfer_strategy(update_child_position);
472
473 ops.emplace_back(collapse);
474 ops_name.emplace_back("COLLAPSE");
475
477 // 1) EdgeSwap
479
480 // swap56
481
482 auto swap56 = std::make_shared<MinOperationSequence>(position_mesh);
483 for (int i = 0; i < 5; ++i) {
484 auto swap = std::make_shared<TetEdgeSwap>(position_mesh, i);
485 swap->collapse().add_invariant(invariant_separate_substructures);
486 swap->collapse().add_invariant(link_condition);
487 swap->collapse().set_new_attribute_strategy(
488 position_handle,
489 CollapseBasicStrategy::CopyOther);
490 swap->collapse().set_new_attribute_strategy(
491 sizing_field_scalar_handle,
492 CollapseBasicStrategy::CopyOther);
493 swap->split().set_new_attribute_strategy(position_handle);
494 swap->split().set_new_attribute_strategy(sizing_field_scalar_handle);
495
496 swap->split().set_new_attribute_strategy(
497 visited_edge_flag_handle,
498 wmtk::operations::SplitBasicStrategy::None,
499 wmtk::operations::SplitRibBasicStrategy::None);
500 swap->collapse().set_new_attribute_strategy(
501 visited_edge_flag_handle,
502 wmtk::operations::CollapseBasicStrategy::None);
503
504 swap->split().set_new_attribute_strategy(
505 target_edge_length_handle,
506 wmtk::operations::SplitBasicStrategy::Copy,
507 wmtk::operations::SplitRibBasicStrategy::Mean);
508 swap->collapse().set_new_attribute_strategy(
509 target_edge_length_handle,
510 wmtk::operations::CollapseBasicStrategy::None);
511
512 swap->add_invariant(std::make_shared<Swap56EnergyBeforeInvariantDouble>(
513 position_mesh,
514 position_handle.as<double>(),
515 i));
516
517 swap->add_transfer_strategy(amips_update);
518
519 swap->collapse().add_invariant(inversion_invariant);
520
521 swap->collapse().add_invariant(envelope_invariant);
522
523 for (const auto& attr : pass_through_attributes) {
524 swap->split().set_new_attribute_strategy(
525 attr,
526 wmtk::operations::SplitBasicStrategy::None,
527 wmtk::operations::SplitRibBasicStrategy::None);
528 swap->collapse().set_new_attribute_strategy(
529 attr,
530 wmtk::operations::CollapseBasicStrategy::None);
531 }
532
533 swap56->add_operation(swap);
534 }
535
536 auto swap56_energy_check = [&](int64_t idx, const simplex::Simplex& t) -> double {
537 constexpr static PrimitiveType PV = PrimitiveType::Vertex;
538 constexpr static PrimitiveType PE = PrimitiveType::Edge;
539 constexpr static PrimitiveType PF = PrimitiveType::Triangle;
540 constexpr static PrimitiveType PT = PrimitiveType::Tetrahedron;
541
542 auto accessor = position_mesh.create_const_accessor(position_handle.as<double>());
543
544 const Tuple e0 = t.tuple();
545 const Tuple e1 = position_mesh.switch_tuple(e0, PV);
546
547 std::array<Tuple, 5> v;
548 auto iter_tuple = e0;
549 for (int64_t i = 0; i < 5; ++i) {
550 v[i] = position_mesh.switch_tuples(iter_tuple, {PE, PV});
551 iter_tuple = position_mesh.switch_tuples(iter_tuple, {PF, PT});
552 }
553 if (iter_tuple != e0) return 0;
554 assert(iter_tuple == e0);
555
556 // five iterable vertices remap to 0-4 by m_collapse_index, 0: m_collapse_index, 5:
557 // e0, 6: e1
558 std::array<Eigen::Vector3<double>, 7> positions = {
559 {accessor.const_vector_attribute(v[(idx + 0) % 5]),
560 accessor.const_vector_attribute(v[(idx + 1) % 5]),
561 accessor.const_vector_attribute(v[(idx + 2) % 5]),
562 accessor.const_vector_attribute(v[(idx + 3) % 5]),
563 accessor.const_vector_attribute(v[(idx + 4) % 5]),
564 accessor.const_vector_attribute(e0),
565 accessor.const_vector_attribute(e1)}};
566
567 std::array<Eigen::Vector3d, 7> positions_double = {
568 {positions[0],
569 positions[1],
570 positions[2],
571 positions[3],
572 positions[4],
573 positions[5],
574 positions[6]}};
575
576 std::array<std::array<int, 4>, 6> new_tets = {
577 {{{0, 1, 2, 5}},
578 {{0, 2, 3, 5}},
579 {{0, 3, 4, 5}},
580 {{0, 1, 2, 6}},
581 {{0, 2, 3, 6}},
582 {{0, 3, 4, 6}}}};
583
584 double new_energy_max = std::numeric_limits<double>::lowest();
585
586 for (int i = 0; i < 6; ++i) {
587 if (wmtk::utils::wmtk_orient3d(
588 positions[new_tets[i][0]],
589 positions[new_tets[i][1]],
590 positions[new_tets[i][2]],
591 positions[new_tets[i][3]]) > 0) {
592 auto energy = wmtk::function::utils::Tet_AMIPS_energy({{
593 positions_double[new_tets[i][0]][0],
594 positions_double[new_tets[i][0]][1],
595 positions_double[new_tets[i][0]][2],
596 positions_double[new_tets[i][1]][0],
597 positions_double[new_tets[i][1]][1],
598 positions_double[new_tets[i][1]][2],
599 positions_double[new_tets[i][2]][0],
600 positions_double[new_tets[i][2]][1],
601 positions_double[new_tets[i][2]][2],
602 positions_double[new_tets[i][3]][0],
603 positions_double[new_tets[i][3]][1],
604 positions_double[new_tets[i][3]][2],
605 }});
606
607
608 if (energy > new_energy_max) new_energy_max = energy;
609 } else {
610 auto energy = wmtk::function::utils::Tet_AMIPS_energy({{
611 positions_double[new_tets[i][1]][0],
612 positions_double[new_tets[i][1]][1],
613 positions_double[new_tets[i][1]][2],
614 positions_double[new_tets[i][0]][0],
615 positions_double[new_tets[i][0]][1],
616 positions_double[new_tets[i][0]][2],
617 positions_double[new_tets[i][2]][0],
618 positions_double[new_tets[i][2]][1],
619 positions_double[new_tets[i][2]][2],
620 positions_double[new_tets[i][3]][0],
621 positions_double[new_tets[i][3]][1],
622 positions_double[new_tets[i][3]][2],
623 }});
624
625
626 if (energy > new_energy_max) new_energy_max = energy;
627 }
628 }
629
630 return new_energy_max;
631 };
632
633 swap56->set_value_function(swap56_energy_check);
634 swap56->add_invariant(std::make_shared<EdgeValenceInvariant>(position_mesh, 5));
635
636
637 // swap44
638 auto swap44 = std::make_shared<MinOperationSequence>(position_mesh);
639 for (int i = 0; i < 2; ++i) {
640 auto swap = std::make_shared<TetEdgeSwap>(position_mesh, i);
641 swap->collapse().add_invariant(invariant_separate_substructures);
642 swap->collapse().add_invariant(link_condition);
643 swap->collapse().set_new_attribute_strategy(
644 position_handle,
645 CollapseBasicStrategy::CopyOther);
646 swap->collapse().set_new_attribute_strategy(
647 sizing_field_scalar_handle,
648 CollapseBasicStrategy::CopyOther);
649 swap->split().set_new_attribute_strategy(position_handle);
650 swap->split().set_new_attribute_strategy(sizing_field_scalar_handle);
651 swap->split().set_new_attribute_strategy(
652 visited_edge_flag_handle,
653 wmtk::operations::SplitBasicStrategy::None,
654 wmtk::operations::SplitRibBasicStrategy::None);
655 swap->collapse().set_new_attribute_strategy(
656 visited_edge_flag_handle,
657 wmtk::operations::CollapseBasicStrategy::None);
658
659 swap->split().set_new_attribute_strategy(
660 target_edge_length_handle,
661 wmtk::operations::SplitBasicStrategy::Copy,
662 wmtk::operations::SplitRibBasicStrategy::Mean);
663 swap->collapse().set_new_attribute_strategy(
664 target_edge_length_handle,
665 wmtk::operations::CollapseBasicStrategy::None);
666
667 swap->add_transfer_strategy(amips_update);
668
669 swap->add_invariant(std::make_shared<Swap44EnergyBeforeInvariantDouble>(
670 position_mesh,
671 position_handle.as<double>(),
672 i));
673
674 // swap->add_invariant(inversion_invariant);
675 swap->collapse().add_invariant(inversion_invariant);
676
677 swap->collapse().add_invariant(envelope_invariant);
678
679 for (const auto& attr : pass_through_attributes) {
680 swap->split().set_new_attribute_strategy(
681 attr,
682 wmtk::operations::SplitBasicStrategy::None,
683 wmtk::operations::SplitRibBasicStrategy::None);
684 swap->collapse().set_new_attribute_strategy(
685 attr,
686 wmtk::operations::CollapseBasicStrategy::None);
687 }
688
689 swap44->add_operation(swap);
690 }
691
692 auto swap44_energy_check = [&](int64_t idx, const simplex::Simplex& t) -> double {
693 constexpr static PrimitiveType PV = PrimitiveType::Vertex;
694 constexpr static PrimitiveType PE = PrimitiveType::Edge;
695 constexpr static PrimitiveType PF = PrimitiveType::Triangle;
696 constexpr static PrimitiveType PT = PrimitiveType::Tetrahedron;
697
698 auto accessor = position_mesh.create_const_accessor(position_handle.as<double>());
699
700 // get the coords of the vertices
701 // input edge end points
702 const Tuple e0 = t.tuple();
703 const Tuple e1 = position_mesh.switch_tuple(e0, PV);
704 // other four vertices
705 std::array<Tuple, 4> v;
706 auto iter_tuple = e0;
707 for (int64_t i = 0; i < 4; ++i) {
708 v[i] = position_mesh.switch_tuples(iter_tuple, {PE, PV});
709 iter_tuple = position_mesh.switch_tuples(iter_tuple, {PF, PT});
710 }
711
712 if (iter_tuple != e0) return 0;
713 assert(iter_tuple == e0);
714
715 std::array<Eigen::Vector3<double>, 6> positions = {
716 {accessor.const_vector_attribute(v[(idx + 0) % 4]),
717 accessor.const_vector_attribute(v[(idx + 1) % 4]),
718 accessor.const_vector_attribute(v[(idx + 2) % 4]),
719 accessor.const_vector_attribute(v[(idx + 3) % 4]),
720 accessor.const_vector_attribute(e0),
721 accessor.const_vector_attribute(e1)}};
722 std::array<Eigen::Vector3d, 6> positions_double = {
723 {positions[0], positions[1], positions[2], positions[3], positions[4], positions[5]}};
724
725 std::array<std::array<int, 4>, 4> new_tets = {
726 {{{0, 1, 2, 4}}, {{0, 2, 3, 4}}, {{0, 1, 2, 5}}, {{0, 2, 3, 5}}}};
727
728 double new_energy_max = std::numeric_limits<double>::lowest();
729
730 for (int i = 0; i < 4; ++i) {
731 if (wmtk::utils::wmtk_orient3d(
732 positions[new_tets[i][0]],
733 positions[new_tets[i][1]],
734 positions[new_tets[i][2]],
735 positions[new_tets[i][3]]) > 0) {
736 auto energy = wmtk::function::utils::Tet_AMIPS_energy({{
737 positions_double[new_tets[i][0]][0],
738 positions_double[new_tets[i][0]][1],
739 positions_double[new_tets[i][0]][2],
740 positions_double[new_tets[i][1]][0],
741 positions_double[new_tets[i][1]][1],
742 positions_double[new_tets[i][1]][2],
743 positions_double[new_tets[i][2]][0],
744 positions_double[new_tets[i][2]][1],
745 positions_double[new_tets[i][2]][2],
746 positions_double[new_tets[i][3]][0],
747 positions_double[new_tets[i][3]][1],
748 positions_double[new_tets[i][3]][2],
749 }});
750
751 if (energy > new_energy_max) new_energy_max = energy;
752 } else {
753 auto energy = wmtk::function::utils::Tet_AMIPS_energy({{
754 positions_double[new_tets[i][1]][0],
755 positions_double[new_tets[i][1]][1],
756 positions_double[new_tets[i][1]][2],
757 positions_double[new_tets[i][0]][0],
758 positions_double[new_tets[i][0]][1],
759 positions_double[new_tets[i][0]][2],
760 positions_double[new_tets[i][2]][0],
761 positions_double[new_tets[i][2]][1],
762 positions_double[new_tets[i][2]][2],
763 positions_double[new_tets[i][3]][0],
764 positions_double[new_tets[i][3]][1],
765 positions_double[new_tets[i][3]][2],
766 }});
767
768 if (energy > new_energy_max) new_energy_max = energy;
769 }
770 }
771
772 return new_energy_max;
773 };
774
775 swap44->set_value_function(swap44_energy_check);
776 swap44->add_invariant(std::make_shared<EdgeValenceInvariant>(position_mesh, 4));
777
778 // swap 32
779 auto swap32 = std::make_shared<TetEdgeSwap>(position_mesh, 0);
780 swap32->add_invariant(std::make_shared<EdgeValenceInvariant>(position_mesh, 3));
781 swap32->add_invariant(std::make_shared<Swap32EnergyBeforeInvariantDouble>(
782 position_mesh,
783 position_handle.as<double>()));
784
785 swap32->collapse().add_invariant(invariant_separate_substructures);
786 swap32->collapse().add_invariant(link_condition);
787 swap32->collapse().set_new_attribute_strategy(
788 position_handle,
789 CollapseBasicStrategy::CopyOther);
790 swap32->collapse().set_new_attribute_strategy(
791 sizing_field_scalar_handle,
792 CollapseBasicStrategy::CopyOther);
793 // swap32->add_invariant(inversion_invariant);
794 swap32->split().set_new_attribute_strategy(position_handle);
795 swap32->split().set_new_attribute_strategy(sizing_field_scalar_handle);
796 swap32->split().set_new_attribute_strategy(
797 visited_edge_flag_handle,
798 wmtk::operations::SplitBasicStrategy::None,
799 wmtk::operations::SplitRibBasicStrategy::None);
800 swap32->collapse().set_new_attribute_strategy(
801 visited_edge_flag_handle,
802 wmtk::operations::CollapseBasicStrategy::None);
803
804 swap32->split().set_new_attribute_strategy(
805 target_edge_length_handle,
806 wmtk::operations::SplitBasicStrategy::Copy,
807 wmtk::operations::SplitRibBasicStrategy::Mean);
808 swap32->collapse().set_new_attribute_strategy(
809 target_edge_length_handle,
810 wmtk::operations::CollapseBasicStrategy::None);
811
812 swap32->add_transfer_strategy(amips_update);
813
814 // hack
815 swap32->collapse().add_invariant(inversion_invariant);
816 swap32->collapse().add_invariant(envelope_invariant);
817
818 for (const auto& attr : pass_through_attributes) {
819 swap32->split().set_new_attribute_strategy(
820 attr,
821 wmtk::operations::SplitBasicStrategy::None,
822 wmtk::operations::SplitRibBasicStrategy::None);
823 swap32->collapse().set_new_attribute_strategy(
824 attr,
825 wmtk::operations::CollapseBasicStrategy::None);
826 }
827
828 auto swap_all = std::make_shared<OrOperationSequence>(position_mesh);
829 swap_all->add_operation(swap32);
830 swap_all->add_operation(swap44);
831 swap_all->add_operation(swap56);
832 swap_all->add_transfer_strategy(tag_update);
833 swap_all->add_transfer_strategy(energy_filter_update);
834 swap_all->add_invariant(
835 std::make_shared<EnergyFilterInvariant>(position_mesh, energy_filter_handle.as<char>()));
836 swap_all->add_invariant(std::make_shared<InteriorEdgeInvariant>(position_mesh));
837 swap_all->add_invariant(std::make_shared<NoChildMeshAttachingInvariant>(position_mesh));
838
839 swap_all->set_priority(long_edges_first);
840
841 swap_all->add_transfer_strategy(update_child_position);
842
843 ops.push_back(swap_all);
844 ops_name.push_back("EDGE SWAP");
845
847 // 4) Smoothing
849
850 using MeshConstrainPair = ProjectOperation::MeshConstrainPair;
851 std::vector<MeshConstrainPair> mesh_constraint_pairs;
852 mesh_constraint_pairs.emplace_back(surface_position_handle, surface_position_handle);
853
854 auto smoothing = std::make_shared<AMIPSOptimizationSmoothingPeriodic>(
855 periodic_mesh,
856 position_mesh,
857 position_handle);
858 smoothing->add_invariant(inversion_invariant);
859 smoothing->add_transfer_strategy(update_child_position);
860 smoothing->use_random_priority() = true;
861
862 smoothing->add_invariant(envelope_invariant);
863 smoothing->add_invariant(
864 std::make_shared<EnergyFilterInvariant>(position_mesh, energy_filter_handle.as<char>()));
865
866 smoothing->add_transfer_strategy(amips_update);
867 smoothing->add_transfer_strategy(edge_length_update);
868 smoothing->add_transfer_strategy(update_child_position);
869
870 ops.push_back(smoothing);
871 ops_name.push_back("SMOOTHING");
872
874 // debug code
876
877 // int kkk = 0;
878
879 // for (const auto& t : periodic_mesh.get_all(PrimitiveType::Edge)) {
880 // auto edges = periodic_mesh.map_to_child(
881 // position_mesh,
882 // simplex::Simplex(periodic_mesh, PrimitiveType::Edge, t));
883 // auto cnt = edges.size();
884 // if (cnt == 0) {
885 // wmtk::logger().error(
886 // "Edge {} cannot map to any tet in child mesh",
887 // periodic_mesh.id(t, PrimitiveType::Edge));
888 // } else if (cnt > 2) {
889 // wmtk::logger().error(
890 // "Edge {} can map to {} edges in child mesh, {} {} {}",
891 // periodic_mesh.id(t, PrimitiveType::Edge),
892 // cnt,
893 // position_mesh.id(edges[0].tuple(), PrimitiveType::Edge),
894 // position_mesh.id(edges[1].tuple(), PrimitiveType::Edge),
895 // position_mesh.id(edges[2].tuple(), PrimitiveType::Edge));
896 // } else if (cnt == 2) {
897 // wmtk::logger().error(
898 // "Edge {} can map to {} edges in child mesh, {} {}",
899 // periodic_mesh.id(t, PrimitiveType::Edge),
900 // cnt,
901 // position_mesh.id(edges[0].tuple(), PrimitiveType::Edge),
902 // position_mesh.id(edges[1].tuple(), PrimitiveType::Edge));
903 // }
904
905 // kkk++;
906 // }
907
908 // wmtk::logger().info("#Periodic Edges: {}",
909 // periodic_mesh.get_all(PrimitiveType::Edge).size()); wmtk::logger().info("#Edges: {}",
910 // position_mesh.get_all(PrimitiveType::Edge).size());
911
912 // write(position_mesh, output_dir, "position_befire_split", "vertices", 0, true);
913
914 // auto mods = (*split)(simplex::Simplex(
915 // periodic_mesh,
916 // PrimitiveType::Edge,
917 // periodic_mesh.tuple_from_id(PrimitiveType::Edge, 251)));
918
919 // kkk = 0;
920
921 // for (const auto& t : periodic_mesh.get_all(PrimitiveType::Edge)) {
922 // auto edges = periodic_mesh.map_to_child(
923 // position_mesh,
924 // simplex::Simplex(periodic_mesh, PrimitiveType::Edge, t));
925 // auto cnt = edges.size();
926 // if (cnt == 0) {
927 // wmtk::logger().error(
928 // "Edge {} cannot map to any tet in child mesh",
929 // periodic_mesh.id(t, PrimitiveType::Edge));
930 // } else if (cnt > 2) {
931 // wmtk::logger().error(
932 // "Edge {} can map to {} edges in child mesh, {} {} {}",
933 // periodic_mesh.id(t, PrimitiveType::Edge),
934 // cnt,
935 // position_mesh.id(edges[0].tuple(), PrimitiveType::Edge),
936 // position_mesh.id(edges[1].tuple(), PrimitiveType::Edge),
937 // position_mesh.id(edges[2].tuple(), PrimitiveType::Edge));
938 // }
939 // kkk++;
940 // }
941
942 // write(position_mesh, output_dir, "psition_after_split", "vertices", 0, true);
943
944 // exit(0);
945
946
948 // Scheduler
950 write(position_mesh, output_dir, output + "_initial", "vertices", 0, intermediate_output);
951
952 Scheduler scheduler;
953
954 int success = 0;
955
956 double old_max_energy = max_amips;
957 double old_avg_energy = 0;
958
959 for (int64_t i = 0; i < passes; ++i) {
960 logger().info("--------------------------- Pass {} ---------------------------", i);
961
962 SchedulerStats pass_stats;
963 int jj = 0;
964
965 for (auto& op : ops) {
966 for (const auto& t : position_mesh.get_all(position_mesh.top_simplex_type())) {
967 const auto vertices = position_mesh.orient_vertices(t);
968 std::vector<Vector3d> pos;
969 for (int i = 0; i < vertices.size(); ++i) {
970 pos.push_back(position_accessor.const_vector_attribute(vertices[i]));
971 }
972 if (wmtk::utils::wmtk_orient3d(pos[0], pos[1], pos[2], pos[3]) <= 0) {
973 wmtk::logger().error("Flipped tet!");
974 }
975 }
976 logger().info("Executing {} ...", ops_name[jj]);
977 SchedulerStats stats;
978 if (op->primitive_type() == PrimitiveType::Edge) {
979 stats = scheduler.run_operation_on_all(*op, visited_edge_flag_handle.as<char>());
980 } else {
981 stats = scheduler.run_operation_on_all(*op);
982 }
983 pass_stats += stats;
984 logger().info(
985 "Executed {}, {} ops (S/F) {}/{}. Time: collecting: {}, sorting: {}, "
986 "executing: {}",
987 ops_name[jj],
988 stats.number_of_performed_operations(),
989 stats.number_of_successful_operations(),
990 stats.number_of_failed_operations(),
991 stats.collecting_time,
992 stats.sorting_time,
993 stats.executing_time);
994
995 // // debug code
996 // auto pos_size = position_mesh.get_all(PrimitiveType::Tetrahedron).size();
997 // auto peri_size = periodic_mesh.get_all(PrimitiveType::Tetrahedron).size();
998
999 // if (pos_size != peri_size) {
1000 // wmtk::logger().error(
1001 // "Unmatch tet count, periodic: {}, position: {}",
1002 // peri_size,
1003 // pos_size);
1004 // }
1005
1006 // int kkkk = 0;
1007 // for (const auto& t : periodic_mesh.get_all(PrimitiveType::Tetrahedron)) {
1008 // auto cnt = periodic_mesh
1009 // .map_to_child(
1010 // position_mesh,
1011 // simplex::Simplex(periodic_mesh, PrimitiveType::Tetrahedron,
1012 // t))
1013 // .size();
1014 // if (cnt == 0) {
1015 // wmtk::logger().error("Tet {} cannot map to any tet in child mesh", kkkk);
1016 // } else if (cnt > 1) {
1017 // wmtk::logger().error("Tet {} can map to {}} tets in child mesh", kkkk, cnt);
1018 // }
1019 // kkkk++;
1020 // }
1021
1022 // int kkk = 0;
1023
1024 // for (const auto& t : periodic_mesh.get_all(PrimitiveType::Edge)) {
1025 // auto edges = periodic_mesh.map_to_child(
1026 // position_mesh,
1027 // simplex::Simplex(periodic_mesh, PrimitiveType::Edge, t));
1028 // auto cnt = edges.size();
1029 // if (cnt == 0) {
1030 // wmtk::logger().error(
1031 // "Edge {} cannot map to any tet in child mesh",
1032 // periodic_mesh.id(t, PrimitiveType::Edge));
1033 // } else if (cnt > 2) {
1034 // wmtk::logger().error(
1035 // "Edge {} can map to {} edges in child mesh, {} {} {}",
1036 // periodic_mesh.id(t, PrimitiveType::Edge),
1037 // cnt,
1038 // position_mesh.id(edges[0].tuple(), PrimitiveType::Edge),
1039 // position_mesh.id(edges[1].tuple(), PrimitiveType::Edge),
1040 // position_mesh.id(edges[2].tuple(), PrimitiveType::Edge));
1041 // // } else if (cnt == 2) {
1042 // // wmtk::logger().info(
1043 // // "Edge {} can map to {} edges in child mesh, {} {}",
1044 // // periodic_mesh.id(t, PrimitiveType::Edge),
1045 // // cnt,
1046 // // position_mesh.id(edges[0].tuple(), PrimitiveType::Edge),
1047 // // position_mesh.id(edges[1].tuple(), PrimitiveType::Edge));
1048 // }
1049
1050 // kkk++;
1051 // }
1052
1053 // write(position_mesh, output_dir, "psition_after_split", "vertices", 0, true);
1054
1055 // if (i > 1) exit(0);
1056
1057
1058 // debug code end
1059
1060
1061 success = stats.number_of_successful_operations();
1062
1063 // compute energy
1064 double avg_energy = 0;
1065 double max_energy = std::numeric_limits<double>::lowest();
1066 double min_energy = std::numeric_limits<double>::max();
1067 for (const auto& t : position_mesh.get_all(position_mesh.top_simplex_type())) {
1068 double e = amips_accessor.scalar_attribute(t);
1069 max_energy = std::max(max_energy, e);
1070 min_energy = std::min(min_energy, e);
1071 avg_energy += e;
1072 }
1073
1074 avg_energy =
1075 avg_energy / position_mesh.get_all(position_mesh.top_simplex_type()).size();
1076
1077 logger().info(
1078 "Max AMIPS Energy: {}, Min AMIPS Energy: {}, Avg AMIPS Energy: {}",
1079 max_energy,
1080 min_energy,
1081 avg_energy);
1082
1083 ++jj;
1084 }
1085
1086 logger().info(
1087 "Executed {} ops (S/F) {}/{}. Time: collecting: {}, sorting: {}, executing: {}",
1088 pass_stats.number_of_performed_operations(),
1089 pass_stats.number_of_successful_operations(),
1090 pass_stats.number_of_failed_operations(),
1091 pass_stats.collecting_time,
1092 pass_stats.sorting_time,
1093 pass_stats.executing_time);
1094
1095 multimesh::consolidate(periodic_mesh);
1096
1097 write(position_mesh, output_dir, output, "vertices", i + 1, intermediate_output);
1098
1099 double max_energy = std::numeric_limits<double>::lowest();
1100 double min_energy = std::numeric_limits<double>::max();
1101 double avg_energy = 0;
1102 for (const auto& t : position_mesh.get_all(position_mesh.top_simplex_type())) {
1103 double e = amips_accessor.scalar_attribute(t);
1104 max_energy = std::max(max_energy, e);
1105 min_energy = std::min(min_energy, e);
1106 avg_energy += e;
1107 }
1108
1109 avg_energy = avg_energy / position_mesh.get_all(position_mesh.top_simplex_type()).size();
1110
1111 logger().info(
1112 "Max AMIPS Energy: {}, Min AMIPS Energy: {}, Avg AMIPS Energy: {}",
1113 max_energy,
1114 min_energy,
1115 avg_energy);
1116
1117 if (i > 0 && old_max_energy - max_energy < 5e-1 &&
1118 (old_avg_energy - avg_energy) / avg_energy < 0.1) {
1119 wmtk::logger().info("adjusting sizing field ...");
1120
1121 adjust_sizing_field(
1122 position_mesh,
1123 position_handle.as<double>(),
1124 edge_length_handle.as<double>(),
1125 sizing_field_scalar_handle.as<double>(),
1126 amips_handle.as<double>(),
1127 target_edge_length_handle.as<double>(),
1128 visited_vertex_flag_handle.as<char>(),
1129 max_amips_energy,
1130 max_energy,
1131 target_edge_length,
1132 0.00001);
1133
1134 wmtk::logger().info("adjusting sizing field finished");
1135 for (const auto& v : position_mesh.get_all(PrimitiveType::Vertex)) {
1136 energy_filter_accessor.scalar_attribute(v) = char(1);
1137 }
1138 wmtk::logger().info("reset energy filter");
1139 } else {
1140 // wmtk::logger().info("setting energy filter ...");
1141 // set_operation_energy_filter(
1142 // position_mesh,
1143 // position_handle.as<double>(),
1144 // amips_handle.as<double>(),
1145 // energy_filter_handle.as<char>(),
1146 // visited_vertex_flag_handle.as<char>(),
1147 // max_amips_energy,
1148 // max_energy,
1149 // target_edge_length);
1150 // wmtk::logger().info("setting energy filter finished");
1151 }
1152
1153 old_max_energy = max_energy;
1154 old_avg_energy = avg_energy;
1155
1156 // stop at good quality
1157 // if (max_energy <= max_amips_energy) break;
1158 }
1159}
1160
1161void write(
1162 const Mesh& mesh,
1163 const std::string& out_dir,
1164 const std::string& name,
1165 const std::string& vname,
1166 const int64_t index,
1167 const bool intermediate_output)
1168{
1169 if (intermediate_output) {
1170 // write tetmesh
1171 const std::filesystem::path data_dir = "";
1172 wmtk::io::ParaviewWriter writer(
1173 data_dir / (name + "_" + std::to_string(index)),
1174 "vertices",
1175 mesh,
1176 true,
1177 true,
1178 true,
1179 true);
1180 mesh.serialize(writer);
1181 }
1182
1183} // namespace
1184
1185void adjust_sizing_field(
1186 Mesh& m,
1187 const TypedAttributeHandle<double>& coordinate_handle,
1188 const TypedAttributeHandle<double>& edge_length_handle,
1189 const TypedAttributeHandle<double>& sizing_field_scalar_handle,
1190 const TypedAttributeHandle<double>& energy_handle,
1191 const TypedAttributeHandle<double>& target_edge_length_handle,
1192 const TypedAttributeHandle<char>& visited_handle,
1193 const double stop_energy,
1194 const double current_max_energy,
1195 const double initial_target_edge_length,
1196 const double min_target_edge_length)
1197{
1198 if (m.top_simplex_type() != PrimitiveType::Tetrahedron) return;
1199
1200 std::cout << "in here" << std::endl;
1201
1202 const auto coordinate_accessor = m.create_const_accessor<double>(coordinate_handle);
1203 const auto edge_length_accessor = m.create_const_accessor<double>(edge_length_handle);
1204 const auto energy_accessor = m.create_const_accessor<double>(energy_handle);
1205
1206 auto sizing_field_scalar_accessor = m.create_accessor<double>(sizing_field_scalar_handle);
1207 auto target_edge_length_accessor = m.create_accessor<double>(target_edge_length_handle);
1208 auto visited_accessor = m.create_accessor<char>(visited_handle);
1209
1210 const double stop_filter_energy = stop_energy * 0.8;
1211 double filter_energy = std::max(current_max_energy / 100, stop_filter_energy);
1212 filter_energy = std::min(filter_energy, 100.0);
1213
1214 const double recover_scalar = 1.5;
1215 const double refine_scalar = 0.5;
1216 const double min_refine_scalar = min_target_edge_length / initial_target_edge_length;
1217
1218 auto vertices_all = m.get_all(PrimitiveType::Vertex);
1219 int64_t v_cnt = vertices_all.size();
1220
1221 std::vector<Vector3d> centroids;
1222 std::queue<Tuple> v_queue;
1223 // std::vector<double> scale_multipliers(v_cnt, recover_scalar);
1224
1225 // get centroids and initial v_queue
1226 for (const auto& t : m.get_all(PrimitiveType::Tetrahedron)) {
1227 // if (std::cbrt(energy_accessor.const_scalar_attribute(t)) < filter_energy) {
1228 if (energy_accessor.const_scalar_attribute(t) < filter_energy) {
1229 // skip good tets
1230 continue;
1231 }
1232 auto vertices = m.orient_vertices(t);
1233 Vector3d c(0, 0, 0);
1234 for (int i = 0; i < 4; ++i) {
1235 c += coordinate_accessor.const_vector_attribute(vertices[i]);
1236 v_queue.emplace(vertices[i]);
1237 }
1238 centroids.emplace_back(c / 4.0);
1239 }
1240
1241 wmtk::logger().info(
1242 "filter energy: {}, low quality tets num: {}",
1243 filter_energy,
1244 centroids.size());
1245
1246 const double R = initial_target_edge_length * 1.8; // update field radius
1247
1248 for (const auto& v : vertices_all) { // reset visited flag
1249 visited_accessor.scalar_attribute(v) = char(0);
1250 }
1251
1252 // TODO: use efficient data structure
1253 auto get_nearest_dist = [&](const Tuple& v) -> double {
1254 Tuple nearest_tuple;
1255 double min_dist = std::numeric_limits<double>::max();
1256 const Vector3d v_pos = coordinate_accessor.const_vector_attribute(v);
1257 for (const auto& c_pos : centroids) {
1258 double dist = (c_pos - v_pos).norm();
1259 min_dist = std::min(min_dist, dist);
1260 }
1261 return min_dist;
1262 };
1263
1264 while (!v_queue.empty()) {
1265 auto v = v_queue.front();
1266 v_queue.pop();
1267
1268 if (visited_accessor.scalar_attribute(v) == char(1)) continue;
1269 visited_accessor.scalar_attribute(v) = char(1);
1270
1271 double dist = std::max(0., get_nearest_dist(v));
1272
1273 if (dist > R) {
1274 visited_accessor.scalar_attribute(v) = char(0);
1275 continue;
1276 }
1277
1278 double scale_multiplier = std::min(
1279 recover_scalar,
1280 dist / R * (1 - refine_scalar) + refine_scalar); // linear interpolate
1281
1282 auto new_scale = sizing_field_scalar_accessor.scalar_attribute(v) * scale_multiplier;
1283 if (new_scale > 1) {
1284 sizing_field_scalar_accessor.scalar_attribute(v) = 1;
1285 } else if (new_scale < min_refine_scalar) {
1286 sizing_field_scalar_accessor.scalar_attribute(v) = min_refine_scalar;
1287 } else {
1288 sizing_field_scalar_accessor.scalar_attribute(v) = new_scale;
1289 }
1290
1291 // push one ring vertices into the queue
1292 for (const auto& v_one_ring : simplex::link(m, simplex::Simplex::vertex(m, v))
1293 .simplex_vector(PrimitiveType::Vertex)) {
1294 if (visited_accessor.scalar_attribute(v_one_ring) == char(1)) continue;
1295 v_queue.push(v_one_ring.tuple());
1296 }
1297 }
1298
1299 // update the rest
1300 for (const auto& v : vertices_all) {
1301 if (visited_accessor.scalar_attribute(v) == char(1)) continue;
1302 auto new_scale = sizing_field_scalar_accessor.scalar_attribute(v) * 1.5;
1303 if (new_scale > 1) {
1304 sizing_field_scalar_accessor.scalar_attribute(v) = 1;
1305 } else if (new_scale < min_refine_scalar) {
1306 sizing_field_scalar_accessor.scalar_attribute(v) = min_refine_scalar;
1307 } else {
1308 sizing_field_scalar_accessor.scalar_attribute(v) = new_scale;
1309 }
1310 }
1311
1312 // update target edge length
1313 for (const auto& e : m.get_all(PrimitiveType::Edge)) {
1314 target_edge_length_accessor.scalar_attribute(e) =
1315 initial_target_edge_length *
1316 (sizing_field_scalar_accessor.scalar_attribute(e) +
1317 sizing_field_scalar_accessor.scalar_attribute(
1318 m.switch_tuple(e, PrimitiveType::Vertex))) /
1319 2;
1320 }
1321}
1322
1323void set_operation_energy_filter(
1324 Mesh& m,
1325 const TypedAttributeHandle<double>& coordinate_handle,
1326 const TypedAttributeHandle<double>& energy_handle,
1327 const TypedAttributeHandle<char>& energy_filter_handle,
1328 const TypedAttributeHandle<char>& visited_handle,
1329 const double stop_energy,
1330 const double current_max_energy,
1331 const double initial_target_edge_length)
1332{
1333 // two ring version
1334 if (m.top_simplex_type() != PrimitiveType::Tetrahedron) return;
1335
1336 const auto coordinate_accessor = m.create_const_accessor<double>(coordinate_handle);
1337 const auto energy_accessor = m.create_const_accessor<double>(energy_handle);
1338
1339 auto energy_filter_accessor = m.create_accessor<char>(energy_filter_handle);
1340 auto visited_accessor = m.create_accessor<char>(visited_handle);
1341
1342 const double stop_filter_energy = stop_energy * 0.8;
1343 double filter_energy = std::max(current_max_energy / 100, stop_filter_energy);
1344 filter_energy = std::min(filter_energy, 100.0);
1345
1346 auto vertices_all = m.get_all(PrimitiveType::Vertex);
1347
1348 for (const auto& v : vertices_all) { // reset visited flag
1349 visited_accessor.scalar_attribute(v) = char(0);
1350 energy_filter_accessor.scalar_attribute(v) = char(0);
1351 }
1352
1353 // get centroids and initial v_queue
1354 for (const auto& t : m.get_all(PrimitiveType::Tetrahedron)) {
1355 if (energy_accessor.const_scalar_attribute(t) < filter_energy) {
1356 // skip good tets
1357 continue;
1358 }
1359 auto vertices = m.orient_vertices(t);
1360 for (const auto& v : vertices) {
1361 energy_filter_accessor.scalar_attribute(v) = char(1);
1362 for (const auto& vv : simplex::k_ring(m, simplex::Simplex::vertex(m, v), 1)
1363 .simplex_vector(PrimitiveType::Vertex)) {
1364 energy_filter_accessor.scalar_attribute(vv) = char(1);
1365 }
1366 }
1367 }
1368}
1369
1370
1371} // namespace wmtk::components::internal
wmtk::Mesh::register_attribute
attribute::MeshAttributeHandle register_attribute(const std::string &name, PrimitiveType type, int64_t size, bool replace=false, T default_value=T(0))
Definition Mesh_attributes.cpp:11
wmtk::Mesh::serialize
void serialize(MeshWriter &writer, const Mesh *local_root=nullptr) const
Definition Mesh.cpp:93
wmtk::Mesh::get_attribute_handle
attribute::MeshAttributeHandle get_attribute_handle(const std::string &name, const PrimitiveType ptype) const
Definition Mesh.hpp:903
wmtk::Mesh::create_const_accessor
const attribute::Accessor< T, Mesh, D > create_const_accessor(const attribute::MeshAttributeHandle &handle) const
wmtk::Mesh::create_accessor
attribute::Accessor< T, Mesh, D > create_accessor(const attribute::MeshAttributeHandle &handle)
wmtk::Mesh::get_all
std::vector< Tuple > get_all(PrimitiveType type) const
Generate a vector of Tuples from global vertex/edge/triangle/tetrahedron index.
Definition Mesh.cpp:18
wmtk::Mesh::switch_tuples
Tuple switch_tuples(const Tuple &tuple, const ContainerType &op_sequence) const
Definition Mesh.hpp:953
wmtk::Mesh::switch_tuple
virtual Tuple switch_tuple(const Tuple &tuple, PrimitiveType type) const =0
switch the orientation of the Tuple of the given dimension
wmtk::Mesh::orient_vertices
virtual std::vector< Tuple > orient_vertices(const Tuple &t) const =0
wmtk::Mesh::top_simplex_type
PrimitiveType top_simplex_type() const
Definition Mesh.hpp:982
wmtk::Scheduler::run_operation_on_all
SchedulerStats run_operation_on_all(operations::Operation &op)
Definition Scheduler.cpp:33
wmtk::SchedulerStats::number_of_failed_operations
int64_t number_of_failed_operations() const
Returns the number of failed operations performed by the scheduler.
Definition Scheduler.hpp:23
wmtk::SchedulerStats::number_of_performed_operations
int64_t number_of_performed_operations() const
Returns the number of performed operations performed by the scheduler.
Definition Scheduler.hpp:30
wmtk::SchedulerStats::number_of_successful_operations
int64_t number_of_successful_operations() const
Returns the number of successful operations performed by the scheduler.
Definition Scheduler.hpp:16
wmtk::Tuple
The Tuple is the basic navigation tool in our mesh data structure.
Definition Tuple.hpp:19
wmtk::attribute::TypedAttributeHandle
Handle that represents attributes for some mesh.
Definition TypedAttributeHandle.hpp:28
wmtk::operations::composite::ProjectOperation::MeshConstrainPair
std::pair< attribute::MeshAttributeHandle, attribute::MeshAttributeHandle > MeshConstrainPair
Definition ProjectOperation.hpp:16
wmtk::simplex::SimplexCollection::simplex_vector
const std::vector< Simplex > & simplex_vector() const
Return const reference to the simplex vector.
Definition SimplexCollection.hpp:25
wmtk::simplex::Simplex::vertex
static Simplex vertex(const Mesh &m, const Tuple &t)
Definition Simplex.hpp:56
extract_child_mesh_from_tag.hpp
PT
constexpr wmtk::PrimitiveType PT
Definition faces_single_dimension.cpp:6
PF
constexpr wmtk::PrimitiveType PF
Definition faces_single_dimension.cpp:5
get_attributes.hpp
k_ring.hpp
wmtk::components::internal::set_operation_energy_filter
void set_operation_energy_filter(Mesh &m, const TypedAttributeHandle< double > &coordinate_handle, const TypedAttributeHandle< double > &energy_handle, const TypedAttributeHandle< char > &energy_filter_handle, const TypedAttributeHandle< char > &visited_handle, const double stop_energy, const double current_max_energy, const double initial_target_edge_length)
Definition periodic_optimization.cpp:1323
wmtk::components::internal::periodic_optimization
void periodic_optimization(Mesh &periodic_mesh, Mesh &position_mesh, Mesh &surface_mesh, const double target_edge_length, const double max_amips_energy, const int64_t passes, const double envelope_size, const bool intermediate_output, std::vector< attribute::MeshAttributeHandle > &pass_through_attributes, std::string output_dir, std::string output)
Definition periodic_optimization.cpp:112
wmtk::components::internal::adjust_sizing_field
void adjust_sizing_field(Mesh &m, const TypedAttributeHandle< double > &coordinate_handle, const TypedAttributeHandle< double > &edge_length_handle, const TypedAttributeHandle< double > &sizing_field_scalar_handle, const TypedAttributeHandle< double > &energy_handle, const TypedAttributeHandle< double > &target_edge_length_handle, const TypedAttributeHandle< char > &visited_handle, const double stop_energy, const double current_max_energy, const double initial_target_edge_length, const double min_target_edge_length)
Definition periodic_optimization.cpp:1185
wmtk::components::internal::write
void write(const Mesh &mesh, const std::string &out_dir, const std::string &name, const std::string &vname, const int64_t index, const bool intermediate_output)
Definition periodic_optimization.cpp:1161
wmtk::function::utils::Tet_AMIPS_energy
double Tet_AMIPS_energy(const std::array< double, 12 > &T)
Definition amips.cpp:380
wmtk::function::Tet_AMIPS_energy
double Tet_AMIPS_energy(const std::array< double, 12 > &T)
Definition AMIPS.cpp:375
wmtk::multimesh::consolidate
void consolidate(Mesh &mesh)
Definition consolidate.cpp:13
wmtk::simplex::vertices
std::vector< Tuple > vertices(const Mesh &m, const Simplex &simplex)
Definition faces_single_dimension.cpp:9
wmtk::simplex::link_condition
bool link_condition(const EdgeMesh &mesh, const Tuple &edge)
Check if the edge to collapse satisfying the link condition.
Definition link_condition.cpp:43
wmtk::simplex::link
SimplexCollection link(const Mesh &mesh, const simplex::Simplex &simplex, const bool sort_and_clean)
Definition link.cpp:84
wmtk::simplex::k_ring
SimplexCollection k_ring(const Mesh &mesh, const Simplex &simplex, int64_t k)
Definition k_ring.cpp:6
wmtk::utils::wmtk_orient3d
int wmtk_orient3d(const Eigen::Ref< const Eigen::Vector3< Rational > > &p0, const Eigen::Ref< const Eigen::Vector3< Rational > > &p1, const Eigen::Ref< const Eigen::Vector3< Rational > > &p2, const Eigen::Ref< const Eigen::Vector3< Rational > > &p3)
Definition orient.cpp:75
wmtk::PE
constexpr PrimitiveType PE
Definition EnvelopeInvariant.cpp:25
wmtk::logger
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:58
wmtk::Vector3d
Vector< double, 3 > Vector3d
Definition Types.hpp:39
wmtk::PV
constexpr PrimitiveType PV
Definition EnvelopeInvariant.cpp:24
periodic_optimization.hpp
data_dir
const std::filesystem::path data_dir
Definition test_simplicial_embedding.cpp:14