MAST3
Multidisciplinary-design Adaptation and Sensitivity Toolkit (MAST)
geometric_filter.hpp
Go to the documentation of this file.
1 /*
2  * MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
3  * Copyright (C) 2013-2020 Manav Bhatia and MAST authors
4  *
5  * This library is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU Lesser General Public
7  * License as published by the Free Software Foundation; either
8  * version 2.1 of the License, or (at your option) any later version.
9  *
10  * This library is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13  * Lesser General Public License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public
16  * License along with this library; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18  */
19 
20 
21 #ifndef __mast__libmesh_geometric_filter_h__
22 #define __mast__libmesh_geometric_filter_h__
23 
24 
25 // MAST includes
26 #include <mast/base/mast_data_types.h>
27 #include <mast/base/exceptions.hpp>
28 #include <mast/numerics/utility.hpp>
29 #include <mast/optimization/design_parameter_vector.hpp>
30 #include <mast/mesh/libmesh/geometric_filter_augment_send_list.hpp>
31 #include <mast/mesh/libmesh/utility.hpp>
32 
33 // libMesh includes
34 #include "libmesh/system.h"
35 #include "libmesh/node.h"
36 #include "libmesh/elem.h"
37 #include "libmesh/mesh_base.h"
38 #include "libmesh/petsc_vector.h"
39 #ifdef LIBMESH_HAVE_NANOFLANN
40 #include "libmesh/nanoflann.hpp"
41 #endif
42 
43 
44 
45 namespace MAST {
46 namespace Mesh {
47 namespace libMeshWrapper {
48 
49 
54 class GeometricFilter {
55 
56 public:
57 
62  GeometricFilter(libMesh::System &sys,
63  const real_t radius):
64  _system (sys),
65  _augment_send_list (nullptr) {
66 
67  _init_filter_radius();
68 
69 #ifdef LIBMESH_HAVE_NANOFLANN
70  _init2(); // KD-tree search using NanoFlann
71 #else
72  _init(); // linear filter search
73 #endif
74 
75  // now initialize and attach sendlist to the dofmap
76  _augment_send_list =
77  new MAST::Mesh::libMeshWrapper::GeometricFilterAugmentSendList(_forward_send_list);
78 
79  _system.get_dof_map().attach_extra_send_list_object(*_augment_send_list);
80 
81  // now we tell the function to
82  _system.get_dof_map().reinit_send_list(_system.get_mesh());
83 
84  // now clear the node list to save memory
85  _nodes.clear();
86  _node_id_to_vec_index.clear();
87  }
88 
89 
90  virtual ~GeometricFilter() {
91 
92  if (_augment_send_list) delete _augment_send_list;
93 
94  // delete the matrix
95  MatDestroy(&_weight_matrix);
96  }
97 
98 
99  inline void disable_sendlist() {
100 
101  Assert0(_augment_send_list, "Send list object not initialized");
102 
103  _augment_send_list->if_enable(false);
104  }
105 
106  template <typename Vec1Type,
107  typename Vec2Type>
108  inline void
109  compute_filtered_values
110  (const Vec1Type &input,
111  Vec2Type &output) const {
112 
113  Assert1(_system.comm().size() == 1,
114  _system.comm().size(),
115  "Method only for scalar runs");
116  Assert2(input.size() == _system.n_dofs(),
117  input.size(), _system.n_dofs(),
118  "Incompatible vector sizes");
119  Assert2(output.size() == _system.n_dofs(),
120  output.size(), _system.n_dofs(),
121  "Incompatible vector sizes");
122 
123  MAST::Numerics::Utility::setZero(output);
124 
125  const uint_t
126  first_local_dof = _system.get_dof_map().first_dof(_system.comm().rank()),
127  last_local_dof = _system.get_dof_map().end_dof(_system.comm().rank());
128 
129  std::map<uint_t, std::vector<std::pair<uint_t, real_t>>>::const_iterator
130  map_it = _filter_map.begin(),
131  map_end = _filter_map.end();
132 
133  for ( ; map_it != map_end; map_it++) {
134 
135  // The forward map only processes the local dofs.
136  if (map_it->first >= first_local_dof &&
137  map_it->first < last_local_dof) {
138 
139  std::vector<std::pair<uint_t, real_t>>::const_iterator
140  vec_it = map_it->second.begin(),
141  vec_end = map_it->second.end();
142 
143  for ( ; vec_it != vec_end; vec_it++) {
144  //if (dvs.is_design_parameter_index(map_it->first))
145  MAST::Numerics::Utility::add
146  (output, map_it->first,
147  MAST::Numerics::Utility::get(input, vec_it->first) * vec_it->second);
148  /*else
149  MAST::Numerics::Utility::set
150  (output, map_it->first,
151  MAST::Numerics::Utility::get(input, vec_it->first));*/
152  }
153  }
154  }
155  }
156 
157 
158  inline void
159  compute_filtered_values(Vec input,
160  Vec output) const {
161 
162 // Assert2(input.size() == _system.n_dofs(),
163 // input.size(), _system.n_dofs(),
164 // "Incompatible vector sizes");
165 // Assert2(output.size() == _system.n_dofs(),
166 // output.size(), _system.n_dofs(),
167 // "Incompatible vector sizes");
168 
169  //MAST::Numerics::Utility::setZero(output);
170 
171  MatMult(_weight_matrix, input, output);
172  }
173 
174  inline void
175  compute_filtered_values(libMesh::NumericVector<real_t> &input,
176  libMesh::NumericVector<real_t> &output) const {
177 
178 // Assert2(input.size() == _system.n_dofs(),
179 // input.size(), _system.n_dofs(),
180 // "Incompatible vector sizes");
181 // Assert2(output.size() == _system.n_dofs(),
182 // output.size(), _system.n_dofs(),
183 // "Incompatible vector sizes");
184 
185  //MAST::Numerics::Utility::setZero(output);
186 
187  this->compute_filtered_values
188  (dynamic_cast<libMesh::PetscVector<real_t>&>(input).vec(),
189  dynamic_cast<libMesh::PetscVector<real_t>&>(output).vec());
190  }
191 
192 
197  template <typename Vec1Type,
198  typename Vec2Type>
199  inline void
200  compute_reverse_filtered_values(const Vec1Type &input,
201  Vec2Type &output) const {
202 
203  Assert1(_system.comm().size() == 1,
204  _system.comm().size(),
205  "Method only for scalar runs");
206  Assert2(input.size() == _system.n_dofs(),
207  input.size(), _system.n_dofs(),
208  "Incompatible vector sizes");
209  Assert2(output.size() == _system.n_dofs(),
210  output.size(), _system.n_dofs(),
211  "Incompatible vector sizes");
212 
213  MAST::Numerics::Utility::setZero(output);
214 
215  const libMesh::DofMap
216  &dof_map = _system.get_dof_map();
217 
218  std::map<uint_t, std::vector<std::pair<uint_t, real_t>>>::const_iterator
219  map_it = _reverse_map.begin(),
220  map_end = _reverse_map.end();
221 
222  for ( ; map_it != map_end; map_it++) {
223 
224  // The reverse map requires processing of dofs that are local and in the send
225  // list. Hence, we use dof_map to check if the dof falls in this category.
226  if (dof_map.semilocal_index(map_it->first)) {
227 
228  std::vector<std::pair<uint_t, real_t>>::const_iterator
229  vec_it = map_it->second.begin(),
230  vec_end = map_it->second.end();
231 
232  for ( ; vec_it != vec_end; vec_it++) {
233  //if (dvs.is_design_parameter_index(map_it->first))
234  MAST::Numerics::Utility::add
235  (output, map_it->first,
236  MAST::Numerics::Utility::get(input, vec_it->first) * vec_it->second);
237  /*else
238  MAST::Numerics::Utility::set
239  (output, map_it->first,
240  MAST::Numerics::Utility::get(input, vec_it->first));*/
241  }
242  }
243  }
244  }
245 
246 
251  inline void
252  compute_reverse_filtered_values(Vec input,
253  Vec output) const {
254 
255 // Assert2(input.size() == _system.n_dofs(),
256 // input.size(), _system.n_dofs(),
257 // "Incompatible vector sizes");
258 // Assert2(output.size() == _system.n_dofs(),
259 // output.size(), _system.n_dofs(),
260 // "Incompatible vector sizes");
261 
262  //MAST::Numerics::Utility::setZero(output);
263 
264  MatMultTranspose(_weight_matrix, input, output);
265  }
266 
267 
268  inline void
269  compute_reverse_filtered_values(libMesh::NumericVector<real_t> &input,
270  libMesh::NumericVector<real_t> &output) const {
271 
272 // Assert2(input.size() == _system.n_dofs(),
273 // input.size(), _system.n_dofs(),
274 // "Incompatible vector sizes");
275 // Assert2(output.size() == _system.n_dofs(),
276 // output.size(), _system.n_dofs(),
277 // "Incompatible vector sizes");
278 
279  //MAST::Numerics::Utility::setZero(output);
280 
281  this->compute_reverse_filtered_values
282  (dynamic_cast<libMesh::PetscVector<real_t>&>(input).vec(),
283  dynamic_cast<libMesh::PetscVector<real_t>&>(output).vec());
284  }
285 
286 
287 
288 
289 private:
290 
291 #ifdef LIBMESH_HAVE_NANOFLANN
292  // Nanoflann uses "duck typing" to allow users to define their own adaptors...
293  template <uint_t Dim>
294  class NanoflannMeshAdaptor
295  {
296  private:
297  // Constant reference to the Mesh we are adapting for use in Nanoflann
298  const libMesh::MeshBase & _mesh;
299  uint_t _n_local_nodes;
300  uint_t _n_nodes;
301 
302  public:
303  NanoflannMeshAdaptor (const libMesh::MeshBase &mesh,
304  const std::vector<const libMesh::Node*> &nodes) :
305  _mesh (mesh),
306  _n_local_nodes (nodes.size()),
307  _n_nodes (mesh.n_nodes()),
308  _nodes (nodes)
309  { }
310 
314  typedef real_t coord_t;
315 
319  inline size_t
320  kdtree_get_point_count() const { return _n_local_nodes; }
321 
326  inline coord_t
327  kdtree_distance(const coord_t * p1,
328  const size_t idx_p2,
329  size_t size) const {
330 
331  Assert2(size == Dim, size, Dim, "Incompatible dimension");
332 
333  // Construct a libmesh Point object from the input coord_t. This
334  // assumes LIBMESH_DIM==3.
335  libMesh::Point point1(p1[0],
336  size > 1 ? p1[1] : 0.,
337  size > 2 ? p1[2] : 0.);
338 
339  // Get the referred-to point from the Mesh
340  const libMesh::Point & point2 = _mesh.point(idx_p2);
341 
342  // Compute Euclidean distance
343  return (point1 - point2).norm_sq();
344  }
345 
351  inline coord_t
352  kdtree_get_pt(const size_t idx, int dim) const
353  {
354  Assert2(dim < (int) Dim, dim, (int) Dim,
355  "Incompatible dimension");
356  Assert2(idx < _n_local_nodes, idx, _n_nodes,
357  "Invalid node index");
358  Assert1(dim < 3, dim, "Invalid dimension");
359 
360  return (*_nodes[idx])(dim);
361  }
362 
369  template <class BBOX>
370  inline bool kdtree_get_bbox(BBOX & /* bb */) const { return false; }
371 
372  const std::vector<const libMesh::Node*> &_nodes;
373  };
374 
375  // Declare a type templated on NanoflannMeshAdaptor
376  typedef nanoflann::L2_Simple_Adaptor<real_t, NanoflannMeshAdaptor<3> > adatper_t;
377  // Declare a KDTree type based on NanoflannMeshAdaptor
378  typedef nanoflann::KDTreeSingleIndexAdaptor<adatper_t, NanoflannMeshAdaptor<3>, 3> kd_tree_t;
379 
380  struct PointData {
381  PointData():
382  dof_index (-1),
383  local_mesh_index (-1),
384  rank (-1),
385  point ({0., 0., 0.})
386  { }
387 
388  uint_t dof_index;
389  uint_t local_mesh_index;
390  uint_t rank;
391  std::vector<real_t> point;
392  };
393 
394 
395  inline void _init2() {
396 
397  libMesh::MeshBase& mesh = _system.get_mesh();
398 
399  // Loop over nodes to try and detect duplicates. We use nanoflann
400  // for this, inspired by
401  // https://gist.github.com/jwpeterson/7a36f9f794df67d51126#file-detect_slit-cc-L65
402  // which was inspired by nanoflann example in libMesh source:
403  // contrib/nanoflann/examples/pointcloud_adaptor_example.cpp
404 
405  // Build adaptor and tree objects
406  NanoflannMeshAdaptor<3> mesh_adaptor(mesh, _nodes);
407  kd_tree_t kd_tree(3, mesh_adaptor, nanoflann::KDTreeSingleIndexAdaptorParams(/*max leaf=*/10));
408 
409  // Construct the tree
410  kd_tree.buildIndex();
411 
412  // synchronize with the notes
413  // stores the bounding box for each rank
414  std::vector<real_t>
415  rank_boxes(6, 0.);
416 
417  // limits of x
418  rank_boxes[0] = kd_tree.root_bbox[0].low;
419  rank_boxes[1] = kd_tree.root_bbox[0].high;
420  // limits of y
421  rank_boxes[2] = kd_tree.root_bbox[1].low;
422  rank_boxes[3] = kd_tree.root_bbox[1].high;
423  // limits of z
424  rank_boxes[4] = kd_tree.root_bbox[2].low;
425  rank_boxes[5] = kd_tree.root_bbox[2].high;
426 
427  // this will modify rank_boxes with such that each block if 6 values is
428  // the box limit coordinates from each rank
429  mesh.comm().allgather(rank_boxes, true);
430 
431  Assert2(rank_boxes.size() == mesh.comm().size()*6,
432  rank_boxes.size(), mesh.comm().size()*6,
433  "Invalid vector dimension after allgather()");
434 
435 
436  real_t
437  d_12 = 0.,
438  sum = 0.,
439  nd_radius = 0.;
440 
441  std::set<uint_t> send_list;
442 
443  const libMesh::DofMap
444  &dof_map = _system.get_dof_map();
445 
446  std::vector<PointData>
447  point_data;
448 
449  const uint_t
450  size = mesh.comm().size(),
451  rank = mesh.comm().rank(),
452  first_local_dof = dof_map.first_dof(rank),
453  last_local_dof = dof_map.end_dof(rank),
454  sys_num = _system.number();
455 
456  uint_t
457  dof_1,
458  dof_2;
459 
460  std::map<const libMesh::Node*, real_t>
461  node_sum;
462 
463  std::vector<const libMesh::Node*>::const_iterator
464  node_it = _nodes.begin(),
465  node_end = _nodes.end();
466 
467  std::vector<std::vector<const libMesh::Node*>>
468  remote_node_dependency(size);
469 
470  for (; node_it != node_end; node_it++) {
471 
472  const libMesh::Node* node = *node_it;
473  nd_radius = _radius->el(node->dof_number(sys_num, 0, 0));
474 
475  // check if this node is near the bbox of each rank
476  for (uint_t i=0; i<size; i++) {
477 
478  if (i == rank) continue;
479 
480 
481  if (// within _radius of x-box on rank i
482  (*node)(0) >= rank_boxes[6*i+0]-nd_radius &&
483  (*node)(0) <= rank_boxes[6*i+1]+nd_radius &&
484  // within _radius of y-box on rank i
485  (*node)(1) >= rank_boxes[6*i+2]-nd_radius &&
486  (*node)(1) <= rank_boxes[6*i+3]+nd_radius &&
487  // within _radius of z-box on rank i
488  (*node)(2) >= rank_boxes[6*i+4]-nd_radius &&
489  (*node)(2) <= rank_boxes[6*i+5]+nd_radius) {
490 
491  remote_node_dependency[i].push_back(node);
492  }
493  }
494  }
495 
496 
497  // now communicate with all processors about the filtered nodes
498  std::vector<real_t>
499  coords_send(size),
500  radius_send(size);
501 
502  std::vector<std::vector<real_t>>
503  coords_recv(size),
504  radius_recv(size),
505  filter_data_node_loc_send(size),
506  filter_data_node_loc_recv(size);
507 
508  std::vector<std::vector<uint_t>>
509  // this stores the number of indices that each node in coords_send will depend on
510  filter_data_n_indices_send(size),
511  filter_data_n_indices_recv(size),
512  // this stores the system dof_indices of nodes that the nodes in coord_send will
513  // depend on.
514  filter_data_dof_index_send(size),
515  filter_data_dof_index_recv(size);
516 
517 
518  for (uint_t i=0; i<size; i++) {
519 
520  for (uint_t j=0; j<size; j++) {
521 
522  if ( i != j) {
523 
524  if (i == rank) {
525 
526  // here, we pack the (x,y,z) coordinates of nodes that the
527  // ith processor wants the jth processor to check
528  coords_send.resize(remote_node_dependency[j].size()*3, 0.);
529  radius_send.resize(remote_node_dependency[j].size(), 0.);
530 
531  for (uint_t k=0; k<remote_node_dependency[j].size(); k++) {
532 
533  coords_send[k*3 + 0] = (*remote_node_dependency[j][k])(0);
534  coords_send[k*3 + 1] = (*remote_node_dependency[j][k])(1);
535  coords_send[k*3 + 2] = (*remote_node_dependency[j][k])(2);
536  radius_send[k] =
537  _radius->el(remote_node_dependency[j][k]->dof_number
538  (sys_num,0,0));
539  }
540 
541  // send information from i to j
542  mesh.comm().send(j, coords_send);
543  mesh.comm().send(j, radius_send);
544 
545  // clear the vector, since we dont need it
546  coords_send.clear();
547  radius_send.clear();
548  }
549  else if (j == rank) {
550 
551  coords_recv[i].clear();
552  radius_recv[i].clear();
553  // get information from i
554  mesh.comm().receive(i, coords_recv[i]);
555  mesh.comm().receive(i, radius_recv[i]);
556 
557  Assert1(coords_recv[i].size()%3 == 0,
558  coords_recv[i].size(),
559  "coordinates must be a multiple of 3");
560  }
561  }
562  }
563  }
564 
565  // now process the data that other processors need from us
566  for (uint_t i=0; i<size; i++) {
567 
568  // nothing to be done if i = rank.
569  if (i == rank) continue;
570 
571  // three coordinates per node, so number of nodes that ith processor
572  // asked for is size/3.
573  uint_t
574  n_nodes = coords_recv[i].size()/3;
575 
576  // initialize this to the number of nodes received from the processor
577  filter_data_n_indices_send[i].reserve(coords_recv[i].size());
578  // we use a conservative estimate of 27 nodes connected to each node
579  filter_data_dof_index_send[i].reserve(coords_recv[i].size()*27);
580  filter_data_node_loc_send[i].reserve (coords_recv[i].size()*27*3);
581 
582  // prepare data for rank i
583  for (uint_t j=0; j<n_nodes; j++) {
584 
585  std::vector<std::pair<size_t, real_t>>
586  indices_dists;
587  nanoflann::RadiusResultSet<real_t, size_t>
588  resultSet(radius_recv[i][j]*radius_recv[i][j], indices_dists);
589 
590  kd_tree.findNeighbors(resultSet, &coords_recv[i][j*3], nanoflann::SearchParams());
591 
592  // add the information from this search result to the information for this node
593  filter_data_n_indices_send[i].push_back(indices_dists.size());
594 
595  for (unsigned r=0; r<indices_dists.size(); ++r) {
596 
597  const libMesh::Node* nd = _nodes[indices_dists[r].first];
598 
599  // location of the node
600  filter_data_node_loc_send[i].push_back((*nd)(0));
601  filter_data_node_loc_send[i].push_back((*nd)(1));
602  filter_data_node_loc_send[i].push_back((*nd)(2));
603  // dof-index of the node
604  filter_data_dof_index_send[i].push_back(nd->dof_number(_system.number(), 0, 0));
605  }
606  }
607  }
608 
609 
610  // now that all the information is ready, we will communicate it from
611  // j to i processor
612  for (uint_t i=0; i<size; i++) {
613 
614  for (uint_t j=0; j<size; j++) {
615 
616  if (i != j) {
617 
618  if (j == rank) {
619 
620  // send information from j to i
621  mesh.comm().send(i, filter_data_n_indices_send[i]);
622  mesh.comm().send(i, filter_data_dof_index_send[i]);
623  mesh.comm().send(i, filter_data_node_loc_send[i]);
624 
625  // now clear the data since we dont need it any more
626  filter_data_n_indices_send[i].clear();
627  filter_data_dof_index_send[i].clear();
628  filter_data_node_loc_send[i].clear();
629  }
630  else if (i == rank) {
631 
632  // get information from j
633  mesh.comm().receive(j, filter_data_n_indices_recv[j]);
634  mesh.comm().receive(j, filter_data_dof_index_recv[j]);
635  mesh.comm().receive(j, filter_data_node_loc_recv[j]);
636  }
637  }
638  }
639  }
640 
641  // now, we search the local nodes and add to it the information from remote couplings
642  // For every node in the mesh, search the KDtree and find any
643  // nodes at _radius distance from the current
644  // node being searched... this will be added to the .
645  node_it = _nodes.begin();
646 
647  for (; node_it != node_end; node_it++) {
648 
649  const libMesh::Node* node = *node_it;
650 
651  dof_1 = node->dof_number(_system.number(), 0, 0);
652  nd_radius = _radius->el(dof_1);
653 
654  // only local dofs are processed.
655  if (dof_map.semilocal_index(dof_1)) {
656 
657  real_t query_pt[3] = {(*node)(0), (*node)(1), (*node)(2)};
658 
659  std::vector<std::pair<size_t, real_t>>
660  indices_dists;
661  indices_dists.push_back(std::pair<size_t, real_t>
662  (_node_id_to_vec_index[node->id()], 0.));
663  nanoflann::RadiusResultSet<real_t, size_t>
664  resultSet(nd_radius*nd_radius, indices_dists);
665 
666  kd_tree.findNeighbors(resultSet, query_pt, nanoflann::SearchParams());
667 
668  sum = 0.;
669 
670  for (unsigned r=0; r<indices_dists.size(); ++r) {
671 
672  d_12 = std::sqrt(indices_dists[r].second);
673 
674  // the distance of this node should be less than or equal to the
675  // specified search radius
676  Assert2(d_12 <= nd_radius, d_12, nd_radius,
677  "Node distance must be <= search radius");
678 
679  sum += nd_radius - d_12;
680  dof_2 = _nodes[indices_dists[r].first]->dof_number(_system.number(), 0, 0);
681 
682  _filter_map[dof_1].push_back(std::pair<uint_t, real_t>
683  (dof_2, nd_radius - d_12));
684 
685  // add this dof to the local send list
686  if (dof_2 < first_local_dof ||
687  dof_2 >= last_local_dof)
688  send_list.insert(dof_2);
689  }
690 
691  Assert1(sum > 0., sum, "Weight must be > 0.");
692 
693  // record this for later use
694  node_sum[node] = sum;
695  }
696  }
697 
698 
699  // also, check the information from remote processors and include them
700  // in the map and send list.
701  uint_t
702  idx = 0,
703  n_remote_nodes = 0;
704 
705  for (uint_t i=0; i<size; i++) {
706 
707  idx = 0;
708 
709  // number of nodes for which we asked i^th rank to search
710  for (uint_t j=0; j<remote_node_dependency[i].size(); j++) {
711 
712  const libMesh::Node*
713  node = remote_node_dependency[i][j];
714 
715  dof_1 = node->dof_number(sys_num, 0, 0);
716  nd_radius = _radius->el(dof_1);
717 
718  // check the number of nodes on the remote node that this
719  // node's filtered values will depend on
720  n_remote_nodes = filter_data_n_indices_recv[i][j];
721 
722  // sum from the local nodes for the geometric filter for
723  // the present node
724  sum = node_sum[node];
725 
726  // now, we will add this information to the respective map
727  for (uint_t k=0; k<n_remote_nodes; k++) {
728 
729  // dof id of this remote node
730  dof_2 = filter_data_dof_index_recv[i][idx];
731 
732  // distance between the present and remote node
733  d_12 =
734  sqrt(pow((*node)(0) - filter_data_node_loc_recv[i][idx*3+0], 2)+
735  pow((*node)(1) - filter_data_node_loc_recv[i][idx*3+1], 2)+
736  pow((*node)(2) - filter_data_node_loc_recv[i][idx*3+2], 2));
737 
738  // contribution to the sum
739  sum += nd_radius - d_12;
740 
741  // add information to dof_1 node about this remote node
742  _filter_map[dof_1].push_back(std::pair<uint_t, real_t>
743  (dof_2, nd_radius - d_12));
744 
745  // add this dof to the local send list
746  send_list.insert(dof_2);
747 
748  // increment the counter for the next data
749  idx++;
750  }
751 
752  // update the value of the sum for this node
753  node_sum[node] = sum;
754  }
755  }
756 
757  // now normalize the weights with respect to the sum
758  // with the coefficients computed for dof_1, divide each coefficient
759  // with the sum
760  node_it = _nodes.begin();
761 
762  for (; node_it != node_end; node_it++) {
763 
764  // node for which the normalization is being done
765  const libMesh::Node* node = *node_it;
766 
767  // dof_id for this node
768  dof_1 = node->dof_number(_system.number(), 0, 0);
769 
770  // filter coefficients for this node
771  std::vector<std::pair<uint_t, real_t>>& vec = _filter_map[dof_1];
772 
773  // sum from the local nodes for the geometric filter for
774  // the present node
775  sum = node_sum[node];
776 
777  for (uint_t i=0; i<vec.size(); i++) {
778 
779  vec[i].second /= sum;
780  Assert1(vec[i].second <= 1., vec[i].second,
781  "Normalized weight must be <= 1.");
782  }
783  }
784 
785  // use the prepared filtering data to initialize the sparse matrix.
786  _init_filter_matrix(_filter_map);
787 
788 
789  // now prepare the reverse map. The send list is sorted for later use.
790  std::set<uint_t>::const_iterator
791  s_it = send_list.begin(),
792  s_end = send_list.end();
793 
794  _forward_send_list.reserve(send_list.size());
795  for ( ; s_it != s_end; s_it++) _forward_send_list.push_back(*s_it);
796 
797  _init_reverse_map(_filter_map, _reverse_map);
798 
799  // compute the largest element size
800  libMesh::MeshBase::const_element_iterator
801  e_it = mesh.local_elements_begin(),
802  e_end = mesh.local_elements_end();
803 
804  for ( ; e_it != e_end; e_it++) {
805  const libMesh::Elem* e = *e_it;
806  d_12 = e->hmax();
807 
808  if (_fe_size < d_12)
809  _fe_size = d_12;
810  }
811 
812  _system.comm().min(_fe_size);
813  }
814 #endif
815 
816 
820  inline void _init() {
821 
822  libMesh::MeshBase& mesh = _system.get_mesh();
823 
824  // currently implemented for replicated mesh
825  Assert0(mesh.is_replicated() ||
826  mesh.comm().size(),
827  "Function implemented only for replicated mesh or serial runs");
828 
829 
830  // iterate over all nodes to compute the
831  libMesh::MeshBase::const_node_iterator
832  node_it_1 = mesh.nodes_begin(),
833  node_it_2 = mesh.nodes_begin(),
834  node_end = mesh.nodes_end();
835 
836  libMesh::Point
837  d;
838 
839  real_t
840  d_12 = 0.,
841  sum = 0.,
842  nd_radius = 0.;
843 
844  std::set<uint_t> send_list;
845 
846  const libMesh::DofMap
847  &dof_map = _system.get_dof_map();
848 
849  const uint_t
850  first_local_dof = _system.get_dof_map().first_dof(_system.comm().rank()),
851  last_local_dof = _system.get_dof_map().end_dof(_system.comm().rank());
852 
853  uint_t
854  sys_num = _system.number(),
855  dof_1,
856  dof_2;
857 
858  for ( ; node_it_1 != node_end; node_it_1++) {
859 
860  dof_1 = (*node_it_1)->dof_number(sys_num, 0, 0);
861  nd_radius = _radius->el(dof_1);
862 
863  // only local dofs are processed.
864  if (dof_map.semilocal_index(dof_1)) {
865 
866  node_it_2 = mesh.nodes_begin();
867  sum = 0.;
868 
869  for ( ; node_it_2 != node_end; node_it_2++) {
870 
871  // compute the distance between the two nodes
872  d = (**node_it_1) - (**node_it_2);
873  d_12 = d.norm();
874 
875  // if the nodes is within the filter radius, add it to the map
876  if (d_12 <= nd_radius) {
877 
878  sum += nd_radius - d_12;
879  dof_2 = (*node_it_2)->dof_number(sys_num, 0, 0);
880 
881  _filter_map[dof_1].push_back(std::pair<uint_t, real_t>(dof_2, nd_radius - d_12));
882 
883  // add this dof to the local send list if it is not a local dof
884  if (dof_2 < first_local_dof ||
885  dof_2 >= last_local_dof)
886  send_list.insert(dof_2);
887  }
888  }
889 
890  Assert1(sum > 0., sum, "Weight must be > 0.");
891 
892  // with the coefficients computed for dof_1, divide each coefficient
893  // with the sum
894  std::vector<std::pair<uint_t, real_t>>& vec = _filter_map[dof_1];
895  for (uint_t i=0; i<vec.size(); i++) {
896 
897  vec[i].second /= sum;
898  Assert1(vec[i].second <= 1., vec[i].second,
899  "Normalized weight must be <= 1.");
900  }
901  }
902  }
903 
904  // use the prepared filtering data to initialize the sparse matrix.
905  _init_filter_matrix(_filter_map);
906 
907 
908  // now prepare the reverse map. The send list is sorted for later use.
909  std::set<uint_t>::const_iterator
910  s_it = send_list.begin(),
911  s_end = send_list.end();
912 
913  _forward_send_list.reserve(send_list.size());
914  for ( ; s_it != s_end; s_it++) _forward_send_list.push_back(*s_it);
915 
916  _init_reverse_map(_filter_map, _reverse_map);
917 
918  // compute the largest element size
919  libMesh::MeshBase::const_element_iterator
920  e_it = mesh.elements_begin(),
921  e_end = mesh.elements_end();
922 
923  for ( ; e_it != e_end; e_it++) {
924  const libMesh::Elem* e = *e_it;
925  d_12 = e->hmax();
926 
927  if (_fe_size < d_12)
928  _fe_size = d_12;
929  }
930  }
931 
932 
933  inline void
934  _init_reverse_map(const std::map<uint_t, std::vector<std::pair<uint_t, real_t>>> &forward_map,
935  std::map<uint_t, std::vector<std::pair<uint_t, real_t>>> &reverse_map) {
936 
937  // now prepare the reverse map
938  std::map<uint_t, std::vector<std::pair<uint_t, real_t>>>::const_iterator
939  it = forward_map.begin(),
940  end = forward_map.end();
941 
942  for ( ; it!=end; it++) {
943 
944  const std::vector<std::pair<uint_t, real_t>>
945  &vec = it->second;
946 
947  for (uint_t i=0; i<vec.size(); i++)
948  reverse_map[vec[i].first].push_back(std::pair<uint_t, real_t>(it->first, vec[i].second));
949  }
950  }
951 
952 
953 
954  inline void
955  _init_filter_radius() {
956 
957  std::unique_ptr<libMesh::NumericVector<real_t>>
958  n_elem_sum(_system.solution->zero_clone().release());
959  _radius.reset(_system.solution->zero_clone().release());
960 
961  // nominal radius is equal to the average element size about a node
962  libMesh::MeshBase::const_element_iterator
963  it = _system.get_mesh().active_local_elements_begin(),
964  end = _system.get_mesh().active_local_elements_end();
965 
966  uint_t
967  first_local_dof = _system.get_dof_map().first_dof(_system.comm().rank()),
968  last_local_dof = _system.get_dof_map().end_dof(_system.comm().rank()),
969  sys_num = _system.number(),
970  dof_num = 0;
971 
972 
973  _nodes.clear();
974  _nodes.reserve(_system.get_mesh().n_local_nodes());
975 
976  std::set<const libMesh::Node*>
977  all_local_nodes;
978 
979  real_t
980  elem_h = 0.;
981 
982  for ( ; it!=end; it++) {
983 
984  const libMesh::Elem *e = *it;
985  elem_h = e->hmax();
986 
987  for (uint_t i=0; i<MAST::Mesh::libMeshWrapper::Utility::n_linear_basis_nodes_on_elem(*e); i++) {
988 
989  const libMesh::Node* nd = e->node_ptr(i);
990  dof_num = nd->dof_number(sys_num, 0, 0);
991 
992  _radius->add(dof_num, elem_h);
993  n_elem_sum->add(dof_num, 1.);
994 
995  // if the node is a local node, add it to the set
996  if (dof_num >= first_local_dof &&
997  dof_num < last_local_dof)
998  all_local_nodes.insert(nd);
999  }
1000  }
1001 
1002  _radius->close();
1003  n_elem_sum->close();
1004 
1005  // now compute the average
1006  for (uint_t i=first_local_dof; i<last_local_dof; i++) {
1007  _radius->set(i, _radius->el(i)/n_elem_sum->el(i));
1008  }
1009  _radius->close();
1010 
1011  // initialize the local nodes
1012  std::set<const libMesh::Node*>::const_iterator
1013  n_it = all_local_nodes.begin(),
1014  n_end = all_local_nodes.end();
1015 
1016  for ( ; n_it != n_end; n_it++) {
1017  _nodes.push_back(*n_it);
1018  _node_id_to_vec_index[(*n_it)->id()] = _nodes.size()-1;
1019  }
1020  }
1021 
1022 
1023 
1024  inline void
1025  _init_filter_matrix
1026  (const std::map<uint_t, std::vector<std::pair<uint_t, real_t>>>& filter_map) {
1027 
1028  const uint_t
1029  n_vals = _system.n_dofs(),
1030  n_local = _system.get_dof_map().n_dofs_on_processor(_system.processor_id()),
1031  first_local_dof = _system.get_dof_map().first_dof(_system.comm().rank()),
1032  last_local_dof = _system.get_dof_map().end_dof(_system.comm().rank());
1033 
1034  PetscErrorCode ierr;
1035 
1036  std::vector<int>
1037  n_nz(n_local, 0),
1038  n_oz(n_local, 0);
1039 
1040  // initialize the data for the sparsity pattern of the matrix
1041  std::map<uint_t, std::vector<std::pair<uint_t, real_t>>>::const_iterator
1042  it = filter_map.begin(),
1043  end = filter_map.end();
1044 
1045  for ( ; it!=end; it++) {
1046 
1047  const std::vector<std::pair<uint_t, real_t>>
1048  &vec = it->second;
1049 
1050  for (uint_t i=0; i<vec.size(); i++) {
1051  uint_t
1052  dof_id = vec[i].first;
1053  if (dof_id >= first_local_dof && dof_id < last_local_dof)
1054  n_nz[it->first - first_local_dof]++;
1055  else
1056  n_oz[it->first - first_local_dof]++;
1057  }
1058  }
1059 
1060 
1061  ierr = MatCreate(_system.comm().get(), &_weight_matrix);
1062  CHKERRABORT(_system.comm().get(), ierr);
1063  ierr = MatSetSizes(_weight_matrix, n_local, n_local, n_vals, n_vals);
1064  CHKERRABORT(_system.comm().get(), ierr);
1065 
1066  std::string nm = _system.name() + "_";
1067  ierr = MatSetOptionsPrefix(_weight_matrix, nm.c_str());
1068  CHKERRABORT(_system.comm().get(), ierr);
1069 
1070  ierr = MatSetFromOptions(_weight_matrix);
1071  CHKERRABORT(_system.comm().get(), ierr);
1072 
1073  ierr = MatSeqAIJSetPreallocation(_weight_matrix,
1074  n_vals,
1075  (PetscInt*)&n_nz[0]);
1076  CHKERRABORT(_system.comm().get(), ierr);
1077  ierr = MatMPIAIJSetPreallocation(_weight_matrix,
1078  0,
1079  (PetscInt*)&n_nz[0],
1080  0,
1081  (PetscInt*)&n_oz[0]);
1082  CHKERRABORT(_system.comm().get(), ierr);
1083  ierr = MatSetOption(_weight_matrix,
1084  MAT_NEW_NONZERO_ALLOCATION_ERR,
1085  PETSC_TRUE);
1086  CHKERRABORT(_system.comm().get(), ierr);
1087 
1088  // now populate the matrix
1089  it = filter_map.begin();
1090  end = filter_map.end();
1091 
1092  for ( ; it!=end; it++) {
1093 
1094  const std::vector<std::pair<uint_t, real_t>>
1095  &vec = it->second;
1096 
1097  for (uint_t i=0; i<vec.size(); i++) {
1098  MatSetValue(_weight_matrix,
1099  it->first,
1100  vec[i].first,
1101  vec[i].second,
1102  INSERT_VALUES);
1103  }
1104  }
1105 
1106  MatAssemblyBegin(_weight_matrix, MAT_FINAL_ASSEMBLY);
1107  CHKERRABORT(_system.comm().get(), ierr);
1108  MatAssemblyEnd(_weight_matrix, MAT_FINAL_ASSEMBLY);
1109  CHKERRABORT(_system.comm().get(), ierr);
1110  }
1111 
1112 
1116  libMesh::System& _system;
1117 
1121  std::unique_ptr<libMesh::NumericVector<real_t>> _radius;
1122 
1123 
1127  real_t _fe_size;
1128 
1133  MAST::Mesh::libMeshWrapper::GeometricFilterAugmentSendList *_augment_send_list;
1134 
1139  std::map<uint_t, std::vector<std::pair<uint_t, real_t>>> _filter_map;
1140 
1144  std::map<uint_t, std::vector<std::pair<uint_t, real_t>>> _reverse_map;
1145 
1149  std::vector<uint_t> _forward_send_list;
1150 
1154  std::vector<const libMesh::Node*> _nodes;
1155 
1156 
1160  std::map<uint_t, uint_t> _node_id_to_vec_index;
1161 
1165  Mat _weight_matrix;
1166 };
1167 
1168 
1169 } // namespace libMeshWrapper
1170 } // namespace Mesh
1171 } // namespace MAST
1172 
1173 
1174 #endif // __mast__libmesh_geometric_filter_h__
MAST::Numerics::Utility::get
ScalarType get(const std::vector< ScalarType > &v, uint_t i)
Definition: utility.hpp:232
MAST::Mesh::libMeshWrapper::GeometricFilter::_init
void _init()
initializes the algebraic data structures
Definition: geometric_filter.hpp:820
geometric_filter_augment_send_list.hpp
MAST::Mesh::libMeshWrapper::GeometricFilter::compute_filtered_values
void compute_filtered_values(const Vec1Type &input, Vec2Type &output) const
Definition: geometric_filter.hpp:110
MAST::Mesh::libMeshWrapper::GeometricFilter::_nodes
std::vector< const libMesh::Node * > _nodes
Local nodes that belong to the first-order Lagrange basis on the elements.
Definition: geometric_filter.hpp:1154
MAST::Numerics::Utility::setZero
void setZero(ValType &m)
Definition: utility.hpp:44
MAST::Mesh::libMeshWrapper::GeometricFilter::_fe_size
real_t _fe_size
largest element size in the level set mesh
Definition: geometric_filter.hpp:1127
Assert1
#define Assert1(cond, v1, msg)
Definition: exceptions.hpp:143
MAST::Mesh::libMeshWrapper::GeometricFilter::_init_reverse_map
void _init_reverse_map(const std::map< uint_t, std::vector< std::pair< uint_t, real_t >>> &forward_map, std::map< uint_t, std::vector< std::pair< uint_t, real_t >>> &reverse_map)
Definition: geometric_filter.hpp:934
MAST::Mesh::libMeshWrapper::GeometricFilter::_radius
std::unique_ptr< libMesh::NumericVector< real_t > > _radius
radius of the filter.
Definition: geometric_filter.hpp:1121
MAST::Mesh::libMeshWrapper::GeometricFilter
Creates a geometric filter for the location-based design variables, for example density and level-set...
Definition: geometric_filter.hpp:54
MAST::Mesh::libMeshWrapper::GeometricFilter::compute_reverse_filtered_values
void compute_reverse_filtered_values(libMesh::NumericVector< real_t > &input, libMesh::NumericVector< real_t > &output) const
Definition: geometric_filter.hpp:269
MAST::Mesh::libMeshWrapper::GeometricFilter::_system
libMesh::System & _system
system on which the level set discrete function is defined
Definition: geometric_filter.hpp:1116
MAST::Mesh::libMeshWrapper::GeometricFilter::_augment_send_list
MAST::Mesh::libMeshWrapper::GeometricFilterAugmentSendList * _augment_send_list
pointer to an object that appends the sendlist for dofmap to localize the dofs needed for local compu...
Definition: geometric_filter.hpp:1133
MAST::Mesh::libMeshWrapper::GeometricFilter::GeometricFilter
GeometricFilter(libMesh::System &sys, const real_t radius)
Definition: geometric_filter.hpp:62
MAST::Mesh::libMeshWrapper::GeometricFilter::_filter_map
std::map< uint_t, std::vector< std::pair< uint_t, real_t > > > _filter_map
Algebraic relation between filtered level set values and the design variables .
Definition: geometric_filter.hpp:1139
mast_data_types.h
MAST::Mesh::libMeshWrapper::GeometricFilter::compute_filtered_values
void compute_filtered_values(Vec input, Vec output) const
Definition: geometric_filter.hpp:159
MAST::Mesh::libMeshWrapper::GeometricFilter::compute_reverse_filtered_values
void compute_reverse_filtered_values(const Vec1Type &input, Vec2Type &output) const
Applies the reverse map, which is used for sensitivity analysis by first computing the sensitivty wrt...
Definition: geometric_filter.hpp:200
MAST::Mesh::libMeshWrapper::GeometricFilter::_node_id_to_vec_index
std::map< uint_t, uint_t > _node_id_to_vec_index
id of local nodes for use in kd-tree search
Definition: geometric_filter.hpp:1160
Assert0
#define Assert0(cond, msg)
Definition: exceptions.hpp:134
MAST::Mesh::libMeshWrapper::Utility::n_linear_basis_nodes_on_elem
uint_t n_linear_basis_nodes_on_elem(const libMesh::Elem &e)
identifies number of ndoes on element
Definition: utility.hpp:39
uint_t
unsigned int uint_t
Definition: mast_data_types.h:43
MAST::Numerics::Utility::add
void add(std::vector< ScalarType > &v, uint_t i, ScalarType s)
Definition: utility.hpp:188
Assert2
#define Assert2(cond, v1, v2, msg)
Definition: exceptions.hpp:152
MAST::Mesh::libMeshWrapper::GeometricFilter::compute_reverse_filtered_values
void compute_reverse_filtered_values(Vec input, Vec output) const
Applies the reverse map, which is used for sensitivity analysis by first computing the sensitivty wrt...
Definition: geometric_filter.hpp:252
design_parameter_vector.hpp
MAST::Mesh::libMeshWrapper::GeometricFilter::_weight_matrix
Mat _weight_matrix
Sparse Matrix object that stores the filtering matrix.
Definition: geometric_filter.hpp:1165
MAST::Mesh::libMeshWrapper::GeometricFilter::_init_filter_matrix
void _init_filter_matrix(const std::map< uint_t, std::vector< std::pair< uint_t, real_t >>> &filter_map)
Definition: geometric_filter.hpp:1026
real_t
double real_t
Definition: mast_data_types.h:42
MAST::Mesh::libMeshWrapper::GeometricFilter::_reverse_map
std::map< uint_t, std::vector< std::pair< uint_t, real_t > > > _reverse_map
this map stores the columns of the matrix, which is required for sensitivity analysis ...
Definition: geometric_filter.hpp:1144
MAST::Mesh::libMeshWrapper::GeometricFilter::_forward_send_list
std::vector< uint_t > _forward_send_list
vector of dof ids that the current processor depends on.
Definition: geometric_filter.hpp:1149
MAST::Mesh::libMeshWrapper::GeometricFilter::compute_filtered_values
void compute_filtered_values(libMesh::NumericVector< real_t > &input, libMesh::NumericVector< real_t > &output) const
Definition: geometric_filter.hpp:175