Wildmeshing Toolkit
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 
71 namespace wmtk::components::internal {
72 
73 using namespace simplex;
74 using namespace operations;
75 using namespace operations::tri_mesh;
76 using namespace operations::tet_mesh;
77 using namespace operations::composite;
78 using namespace function;
79 using namespace invariants;
80 
81 void 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 
89 void 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 
102 void 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 
112 void 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 
1161 void 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 
1185 void 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 
1323 void 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::orient_vertices
virtual std::vector< Tuple > orient_vertices(const Tuple &t) const =0
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:92
wmtk::Mesh::get_attribute_handle
attribute::MeshAttributeHandle get_attribute_handle(const std::string &name, const PrimitiveType ptype) const
Definition: Mesh.hpp:917
wmtk::Mesh::create_const_accessor
const attribute::Accessor< T, Mesh, D > create_const_accessor(const attribute::MeshAttributeHandle &handle) const
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:967
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::create_accessor
attribute::Accessor< T, Mesh, D > create_accessor(const attribute::MeshAttributeHandle &handle)
wmtk::Mesh::top_simplex_type
PrimitiveType top_simplex_type() const
Definition: Mesh.hpp:996
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::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::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::components::output::output
void output(const Mesh &mesh, const std::filesystem::path &file, const std::string &position_attr_name)
Write the mesh to file.
Definition: output.cpp:16
wmtk::function::utils::amips
auto amips(const Eigen::MatrixBase< Derived > &B)
Definition: amips.hpp:20
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::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::vertices
std::vector< Tuple > vertices(const Mesh &m, const Simplex &simplex)
Definition: faces_single_dimension.cpp:9
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::PV
constexpr PrimitiveType PV
Definition: EnvelopeInvariant.cpp:24
wmtk::logger
spdlog::logger & logger()
Retrieves the current logger.
Definition: Logger.cpp:58
wmtk::Vector3d
Vector< double, 3 > Vector3d
Definition: Types.hpp:39
wmtk::PE
constexpr PrimitiveType PE
Definition: EnvelopeInvariant.cpp:25
periodic_optimization.hpp
data_dir
const std::filesystem::path data_dir
Definition: test_component_delaunay.cpp:15