mesh_smoother_vsmoother.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 #include "libmesh/libmesh_config.h"
19 #ifdef LIBMESH_ENABLE_VSMOOTHER
20 
21 // Local includes
22 #include "libmesh/mesh_smoother_vsmoother.h"
23 #include "libmesh/mesh_tools.h"
24 #include "libmesh/elem.h"
25 #include "libmesh/unstructured_mesh.h"
26 #include "libmesh/utility.h"
27 
28 // C++ includes
29 #include <time.h> // for clock_t, clock()
30 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
31 #include <cmath>
32 #include <iomanip>
33 #include <limits>
34 
35 namespace libMesh
36 {
37 
38 // Optimization at -O2 or greater seem to break Intel's icc. So if we are
39 // being compiled with icc let's dumb-down the optimizations for this file
40 #ifdef __INTEL_COMPILER
41 # pragma optimize ( "", off )
42 #endif
43 
44 // Member functions for the Variational Smoother
45 VariationalMeshSmoother::VariationalMeshSmoother(UnstructuredMesh & mesh,
46  double theta,
47  unsigned miniter,
48  unsigned maxiter,
49  unsigned miniterBC) :
50  MeshSmoother(mesh),
51  _percent_to_move(1),
52  _dist_norm(0.),
53  _adapt_data(nullptr),
54  _dim(mesh.mesh_dimension()),
55  _miniter(miniter),
56  _maxiter(maxiter),
57  _miniterBC(miniterBC),
58  _metric(UNIFORM),
59  _adaptive_func(NONE),
60  _theta(theta),
61  _generate_data(false),
62  _n_nodes(0),
63  _n_cells(0),
64  _n_hanging_edges(0),
65  _area_of_interest(nullptr)
66 {}
67 
68 
69 
70 
71 VariationalMeshSmoother::VariationalMeshSmoother(UnstructuredMesh & mesh,
72  std::vector<float> * adapt_data,
73  double theta,
74  unsigned miniter,
75  unsigned maxiter,
76  unsigned miniterBC,
77  double percent_to_move) :
78  MeshSmoother(mesh),
79  _percent_to_move(percent_to_move),
80  _dist_norm(0.),
81  _adapt_data(adapt_data),
82  _dim(mesh.mesh_dimension()),
83  _miniter(miniter),
84  _maxiter(maxiter),
85  _miniterBC(miniterBC),
86  _metric(UNIFORM),
87  _adaptive_func(CELL),
88  _theta(theta),
89  _generate_data(false),
90  _n_nodes(0),
91  _n_cells(0),
92  _n_hanging_edges(0),
93  _area_of_interest(nullptr)
94 {}
95 
96 
97 
98 VariationalMeshSmoother::VariationalMeshSmoother(UnstructuredMesh & mesh,
99  const UnstructuredMesh * area_of_interest,
100  std::vector<float> * adapt_data,
101  double theta,
102  unsigned miniter,
103  unsigned maxiter,
104  unsigned miniterBC,
105  double percent_to_move) :
106  MeshSmoother(mesh),
107  _percent_to_move(percent_to_move),
108  _dist_norm(0.),
109  _adapt_data(adapt_data),
110  _dim(mesh.mesh_dimension()),
111  _miniter(miniter),
112  _maxiter(maxiter),
113  _miniterBC(miniterBC),
114  _metric(UNIFORM),
115  _adaptive_func(CELL),
116  _theta(theta),
117  _generate_data(false),
118  _n_nodes(0),
119  _n_cells(0),
120  _n_hanging_edges(0),
121  _area_of_interest(area_of_interest)
122 {}
123 
124 
125 
126 double VariationalMeshSmoother::smooth(unsigned int)
127 {
128  // If the log file is already open, for example on subsequent calls
129  // to smooth() on the same object, we'll just keep writing to it,
130  // otherwise we'll open it...
131  if (!_logfile.is_open())
132  _logfile.open("smoother.out");
133 
134  int
135  me = _metric,
136  gr = _generate_data ? 0 : 1,
137  adp = _adaptive_func,
138  miniter = _miniter,
139  maxiter = _maxiter,
140  miniterBC = _miniterBC;
141 
142  double theta = _theta;
143 
144  // Metric file name
145  std::string metric_filename = "smoother.metric";
146  if (gr == 0 && me > 1)
147  {
148  // grid filename
149  std::string grid_filename = "smoother.grid";
150 
151  // generate metric from initial mesh (me = 2,3)
152  metr_data_gen(grid_filename, metric_filename, me);
153  }
154 
155  // Initialize the _n_nodes and _n_cells member variables
156  this->_n_nodes = _mesh.n_nodes();
157  this->_n_cells = _mesh.n_active_elem();
158 
159  // Initialize the _n_hanging_edges member variable
160  MeshTools::find_hanging_nodes_and_parents(_mesh, _hanging_nodes);
161  this->_n_hanging_edges =
162  cast_int<dof_id_type>(_hanging_nodes.size());
163 
164  std::vector<int>
165  mask(_n_nodes),
166  edges(2*_n_hanging_edges),
167  mcells(_n_cells),
168  hnodes(_n_hanging_edges);
169 
170  Array2D<double> R(_n_nodes, _dim);
171  Array2D<int> cells(_n_cells, 3*_dim + _dim%2);
172  Array3D<double> H(_n_cells, _dim, _dim);
173 
174  // initial grid
175  int vms_err = readgr(R, mask, cells, mcells, edges, hnodes);
176  if (vms_err < 0)
177  {
178  _logfile << "Error reading input mesh file" << std::endl;
179  return _dist_norm;
180  }
181 
182  if (me > 1)
183  vms_err = readmetr(metric_filename, H);
184 
185  if (vms_err < 0)
186  {
187  _logfile << "Error reading metric file" << std::endl;
188  return _dist_norm;
189  }
190 
191  std::vector<int> iter(4);
192  iter[0] = miniter;
193  iter[1] = maxiter;
194  iter[2] = miniterBC;
195 
196  // grid optimization
197  _logfile << "Starting Grid Optimization" << std::endl;
198  clock_t ticks1 = clock();
199  full_smooth(R, mask, cells, mcells, edges, hnodes, theta, iter, me, H, adp, gr);
200  clock_t ticks2 = clock();
201  _logfile << "full_smooth took ("
202  << ticks2
203  << "-"
204  << ticks1
205  << ")/"
206  << CLOCKS_PER_SEC
207  << " = "
208  << static_cast<double>(ticks2-ticks1)/static_cast<double>(CLOCKS_PER_SEC)
209  << " seconds"
210  << std::endl;
211 
212  // save result
213  _logfile << "Saving Result" << std::endl;
214  writegr(R);
215 
216  libmesh_assert_greater (_dist_norm, 0);
217  return _dist_norm;
218 }
219 
220 
221 
222 // save grid
223 int VariationalMeshSmoother::writegr(const Array2D<double> & R)
224 {
225  libMesh::out << "Starting writegr" << std::endl;
226 
227  // Adjust nodal coordinates to new positions
228  {
229  libmesh_assert_equal_to(_dist_norm, 0.);
230  _dist_norm = 0;
231  int i = 0;
232  for (auto & node : _mesh.node_ptr_range())
233  {
234  double total_dist = 0.;
235 
236  // Get a reference to the node
237  Node & node_ref = *node;
238 
239  // For each node set its X Y [Z] coordinates
240  for (unsigned int j=0; j<_dim; j++)
241  {
242  double distance = R[i][j] - node_ref(j);
243 
244  // Save the squares of the distance
245  total_dist += Utility::pow<2>(distance);
246 
247  node_ref(j) += distance * _percent_to_move;
248  }
249 
250  libmesh_assert_greater_equal (total_dist, 0.);
251 
252  // Add the distance this node moved to the global distance
253  _dist_norm += total_dist;
254 
255  i++;
256  }
257 
258  // Relative "error"
259  _dist_norm = std::sqrt(_dist_norm/_mesh.n_nodes());
260  }
261 
262  libMesh::out << "Finished writegr" << std::endl;
263  return 0;
264 }
265 
266 
267 
268 // reading grid from input file
269 int VariationalMeshSmoother::readgr(Array2D<double> & R,
270  std::vector<int> & mask,
271  Array2D<int> & cells,
272  std::vector<int> & mcells,
273  std::vector<int> & edges,
274  std::vector<int> & hnodes)
275 {
276  libMesh::out << "Starting readgr" << std::endl;
277  // add error messages where format can be inconsistent
278 
279  // Create a quickly-searchable list of boundary nodes.
280  std::unordered_set<dof_id_type> boundary_node_ids =
281  MeshTools::find_boundary_nodes (_mesh);
282 
283  // Grab node coordinates and set mask
284  {
285  // Only compute the node to elem map once
286  std::unordered_map<dof_id_type, std::vector<const Elem *>> nodes_to_elem_map;
287  MeshTools::build_nodes_to_elem_map(_mesh, nodes_to_elem_map);
288 
289  int i = 0;
290  for (auto & node : _mesh.node_ptr_range())
291  {
292  // Get a reference to the node
293  Node & node_ref = *node;
294 
295  // For each node grab its X Y [Z] coordinates
296  for (unsigned int j=0; j<_dim; j++)
297  R[i][j] = node_ref(j);
298 
299  // Set the Proper Mask
300  // Internal nodes are 0
301  // Immovable boundary nodes are 1
302  // Movable boundary nodes are 2
303  if (boundary_node_ids.count(i))
304  {
305  // Only look for sliding edge nodes in 2D
306  if (_dim == 2)
307  {
308  // Find all the nodal neighbors... that is the nodes directly connected
309  // to this node through one edge
310  std::vector<const Node *> neighbors;
311  MeshTools::find_nodal_neighbors(_mesh, node_ref, nodes_to_elem_map, neighbors);
312 
313  // Grab the x,y coordinates
314  Real x = node_ref(0);
315  Real y = node_ref(1);
316 
317  // Theta will represent the atan2 angle (meaning with the proper quadrant in mind)
318  // of the neighbor node in a system where the current node is at the origin
319  Real theta = 0;
320  std::vector<Real> thetas;
321 
322  // Calculate the thetas
323  for (const auto & neighbor : neighbors)
324  {
325  // Note that the x and y values of this node are subtracted off
326  // this centers the system around this node
327  theta = atan2((*neighbor)(1)-y, (*neighbor)(0)-x);
328 
329  // Save it for later
330  thetas.push_back(theta);
331  }
332 
333  // Assume the node is immovable... then prove otherwise
334  mask[i] = 1;
335 
336  // Search through neighbor nodes looking for two that form a straight line with this node
337  for (std::size_t a=0; a<thetas.size()-1; a++)
338  {
339  // Only try each pairing once
340  for (std::size_t b=a+1; b<thetas.size(); b++)
341  {
342  // Find if the two neighbor nodes angles are 180 degrees (pi) off of each other (withing a tolerance)
343  // In order to make this a true movable boundary node... the two that forma straight line with
344  // it must also be on the boundary
345  if (boundary_node_ids.count(neighbors[a]->id()) &&
346  boundary_node_ids.count(neighbors[b]->id()) &&
347  (
348  (std::abs(thetas[a] - (thetas[b] + (libMesh::pi))) < .001) ||
349  (std::abs(thetas[a] - (thetas[b] - (libMesh::pi))) < .001)
350  )
351  )
352  {
353  // if ((*(*it))(1) > 0.25 || (*(*it))(0) > .7 || (*(*it))(0) < .01)
354  mask[i] = 2;
355  }
356 
357  }
358  }
359  }
360  else // In 3D set all boundary nodes to be fixed
361  mask[i] = 1;
362  }
363  else
364  mask[i] = 0; // Internal Node
365 
366  // libMesh::out << "Node: " << i << " Mask: " << mask[i] << std::endl;
367  i++;
368  }
369  }
370 
371  // Grab the connectivity
372  // FIXME: Generalize this!
373  {
374  int i = 0;
375  for (const auto & elem : _mesh.active_element_ptr_range())
376  {
377  // Keep track of the number of nodes
378  // there must be 6 for 2D
379  // 10 for 3D
380  // If there are actually less than that -1 must be used
381  // to fill out the rest
382  int num = 0;
383  /*
384  int num_necessary = 0;
385 
386  if (_dim == 2)
387  num_necessary = 6;
388  else
389  num_necessary = 10;
390  */
391 
392  if (_dim == 2)
393  {
394  switch (elem->n_vertices())
395  {
396  // Grab nodes that do exist
397  case 3: // Tri
398  for (unsigned int k=0; k<elem->n_vertices(); k++)
399  cells[i][k] = elem->node_id(k);
400 
401  num = elem->n_vertices();
402  break;
403 
404  case 4: // Quad 4
405  cells[i][0] = elem->node_id(0);
406  cells[i][1] = elem->node_id(1);
407  cells[i][2] = elem->node_id(3); // Note that 2 and 3 are switched!
408  cells[i][3] = elem->node_id(2);
409  num = 4;
410  break;
411 
412  default:
413  libmesh_error_msg("Unknown number of vertices = " << elem->n_vertices());
414  }
415  }
416  else
417  {
418  // Grab nodes that do exist
419  switch (elem->n_vertices())
420  {
421  // Tet 4
422  case 4:
423  for (unsigned int k=0; k<elem->n_vertices(); k++)
424  cells[i][k] = elem->node_id(k);
425  num = elem->n_vertices();
426  break;
427 
428  // Hex 8
429  case 8:
430  cells[i][0] = elem->node_id(0);
431  cells[i][1] = elem->node_id(1);
432  cells[i][2] = elem->node_id(3); // Note that 2 and 3 are switched!
433  cells[i][3] = elem->node_id(2);
434 
435  cells[i][4] = elem->node_id(4);
436  cells[i][5] = elem->node_id(5);
437  cells[i][6] = elem->node_id(7); // Note that 6 and 7 are switched!
438  cells[i][7] = elem->node_id(6);
439  num=8;
440  break;
441 
442  default:
443  libmesh_error_msg("Unknown 3D element with = " << elem->n_vertices() << " vertices.");
444  }
445  }
446 
447  // Fill in the rest with -1
448  for (int j=num; j<static_cast<int>(cells[i].size()); j++)
449  cells[i][j] = -1;
450 
451  // Mask it with 0 to state that this is an active element
452  // FIXME: Could be something other than zero
453  mcells[i] = 0;
454  i++;
455  }
456  }
457 
458  // Grab hanging node connectivity
459  {
460  std::map<dof_id_type, std::vector<dof_id_type>>::iterator
461  it = _hanging_nodes.begin(),
462  end = _hanging_nodes.end();
463 
464  for (int i=0; it!=end; it++)
465  {
466  libMesh::out << "Parent 1: " << (it->second)[0] << std::endl;
467  libMesh::out << "Parent 2: " << (it->second)[1] << std::endl;
468  libMesh::out << "Hanging Node: " << it->first << std::endl << std::endl;
469 
470  // First Parent
471  edges[2*i] = (it->second)[1];
472 
473  // Second Parent
474  edges[2*i+1] = (it->second)[0];
475 
476  // Hanging Node
477  hnodes[i] = it->first;
478 
479  i++;
480  }
481  }
482  libMesh::out << "Finished readgr" << std::endl;
483 
484  return 0;
485 }
486 
487 
488 
489 // Read Metrics
490 int VariationalMeshSmoother::readmetr(std::string name,
491  Array3D<double> & H)
492 {
493  std::ifstream infile(name.c_str());
494  std::string dummy;
495 
496  for (dof_id_type i=0; i<_n_cells; i++)
497  for (unsigned j=0; j<_dim; j++)
498  {
499  for (unsigned k=0; k<_dim; k++)
500  infile >> H[i][j][k];
501 
502  // Read to end of line and discard
503  std::getline(infile, dummy);
504  }
505 
506  return 0;
507 }
508 
509 
510 
511 // Stolen from ErrorVector!
512 float VariationalMeshSmoother::adapt_minimum() const
513 {
514  float min = std::numeric_limits<float>::max();
515 
516  for (std::size_t i=0; i<_adapt_data->size(); i++)
517  {
518  // Only positive (or zero) values in the error vector
519  libmesh_assert_greater_equal ((*_adapt_data)[i], 0.);
520  min = std::min (min, (*_adapt_data)[i]);
521  }
522 
523  // ErrorVectors are for positive values
524  libmesh_assert_greater_equal (min, 0.);
525 
526  return min;
527 }
528 
529 
530 
531 void VariationalMeshSmoother::adjust_adapt_data()
532 {
533  // For convenience
534  const UnstructuredMesh & aoe_mesh = *_area_of_interest;
535  std::vector<float> & adapt_data = *_adapt_data;
536 
537  float min = adapt_minimum();
538 
539  MeshBase::const_element_iterator el = _mesh.elements_begin();
540  const MeshBase::const_element_iterator end_el = _mesh.elements_end();
541 
542  MeshBase::const_element_iterator aoe_el = aoe_mesh.elements_begin();
543  const MeshBase::const_element_iterator aoe_end_el = aoe_mesh.elements_end();
544 
545  // Counter to keep track of which spot in adapt_data we're looking at
546  for (int i=0; el!=end_el; el++)
547  {
548  // Only do this for active elements
549  if (adapt_data[i])
550  {
551  Point centroid = (*el)->centroid();
552  bool in_aoe = false;
553 
554  // See if the elements centroid lies in the aoe mesh
555  // for (aoe_el=aoe_mesh.elements_begin(); aoe_el != aoe_end_el; ++aoe_el)
556  // {
557  if ((*aoe_el)->contains_point(centroid))
558  in_aoe = true;
559  // }
560 
561  // If the element is not in the area of interest... then set
562  // it's adaptivity value to the minimum.
563  if (!in_aoe)
564  adapt_data[i] = min;
565  }
566  i++;
567  }
568 }
569 
570 
571 
572 // Read Adaptivity
573 int VariationalMeshSmoother::read_adp(std::vector<double> & afun)
574 {
575  std::vector<float> & adapt_data = *_adapt_data;
576 
577  if (_area_of_interest)
578  adjust_adapt_data();
579 
580  std::size_t m = adapt_data.size();
581 
582  std::size_t j =0 ;
583  for (std::size_t i=0; i<m; i++)
584  if (adapt_data[i] != 0)
585  {
586  afun[j] = adapt_data[i];
587  j++;
588  }
589 
590  return 0;
591 }
592 
593 
594 
595 double VariationalMeshSmoother::jac3(double x1,
596  double y1,
597  double z1,
598  double x2,
599  double y2,
600  double z2,
601  double x3,
602  double y3,
603  double z3)
604 {
605  return x1*(y2*z3 - y3*z2) + y1*(z2*x3 - z3*x2) + z1*(x2*y3 - x3*y2);
606 }
607 
608 
609 
610 double VariationalMeshSmoother::jac2(double x1,
611  double y1,
612  double x2,
613  double y2)
614 {
615  return x1*y2 - x2*y1;
616 }
617 
618 
619 
620 // BasisA determines matrix H^(-T)Q on one Jacobian matrix
621 int VariationalMeshSmoother::basisA(Array2D<double> & Q,
622  int nvert,
623  const std::vector<double> & K,
624  const Array2D<double> & H,
625  int me)
626 {
627  Array2D<double> U(_dim, nvert);
628 
629  // Some useful constants
630  const double
631  sqrt3 = std::sqrt(3.),
632  sqrt6 = std::sqrt(6.);
633 
634  if (_dim == 2)
635  {
636  if (nvert == 4)
637  {
638  // quad
639  U[0][0] = -(1-K[1]);
640  U[0][1] = (1-K[1]);
641  U[0][2] = -K[1];
642  U[0][3] = K[1];
643 
644  U[1][0] = -(1-K[0]);
645  U[1][1] = -K[0];
646  U[1][2] = (1-K[0]);
647  U[1][3] = K[0];
648  }
649  else if (nvert == 3)
650  {
651  // tri
652  // for right target triangle
653  // U[0][0] = -1.; U[1][0] = -1.;
654  // U[0][1] = 1.; U[1][1] = 0.;
655  // U[0][2] = 0.; U[1][2] = 1.;
656 
657  // for regular triangle
658  U[0][0] = -1.;
659  U[0][1] = 1.;
660  U[0][2] = 0;
661 
662  U[1][0] = -1./sqrt3;
663  U[1][1] = -1./sqrt3;
664  U[1][2] = 2./sqrt3;
665  }
666  else if (nvert == 6)
667  {
668  // curved triangle
669  U[0][0] = -3*(1-K[0]-K[1]*0.5)+(K[0]-0.5*K[1])+K[1];
670  U[0][1] = -(1-K[0]-K[1]*0.5)+(K[0]-0.5*K[1])*3-K[1];
671  U[0][2] = 0;
672  U[0][3] = K[1]*4;
673  U[0][4] = -4*K[1];
674  U[0][5] = 4*(1-K[0]-K[1]*0.5)-4*(K[0]-0.5*K[1]);
675 
676  U[1][0] = (1-K[2]-K[3]*0.5)*(-1.5)+(K[2]-0.5*K[3])*(0.5)+K[3]*(0.5);
677  U[1][1] = (1-K[2]-K[3]*0.5)*(0.5)+(K[2]-0.5*K[3])*(-1.5)+K[3]*(0.5);
678  U[1][2] = (1-K[2]-K[3]*0.5)*(-1)+(K[2]-0.5*K[3])*(-1)+K[3]*(3);
679  U[1][3] = (K[2]-0.5*K[3])*(4.)+K[3]*(-2.);
680  U[1][4] = (1-K[2]-K[3]*0.5)*(4.)+K[3]*(-2.);
681  U[1][5] = (1-K[2]-K[3]*0.5)*(-2.)+(K[2]-0.5*K[3])*(-2.);
682  }
683  else
684  libmesh_error_msg("Invalid nvert = " << nvert);
685  }
686 
687  if (_dim == 3)
688  {
689  if (nvert == 8)
690  {
691  // hex
692  U[0][0] = -(1-K[0])*(1-K[1]);
693  U[0][1] = (1-K[0])*(1-K[1]);
694  U[0][2] = -K[0]*(1-K[1]);
695  U[0][3] = K[0]*(1-K[1]);
696  U[0][4] = -(1-K[0])*K[1];
697  U[0][5] = (1-K[0])*K[1];
698  U[0][6] = -K[0]*K[1];
699  U[0][7] = K[0]*K[1];
700 
701  U[1][0] = -(1-K[2])*(1-K[3]);
702  U[1][1] = -K[2]*(1-K[3]);
703  U[1][2] = (1-K[2])*(1-K[3]);
704  U[1][3] = K[2]*(1-K[3]);
705  U[1][4] = -(1-K[2])*K[3];
706  U[1][5] = -K[2]*K[3];
707  U[1][6] = (1-K[2])*K[3];
708  U[1][7] = K[2]*K[3];
709 
710  U[2][0] = -(1-K[4])*(1-K[5]);
711  U[2][1] = -K[4]*(1-K[5]);
712  U[2][2] = -(1-K[4])*K[5];
713  U[2][3] = -K[4]*K[5];
714  U[2][4] = (1-K[4])*(1-K[5]);
715  U[2][5] = K[4]*(1-K[5]);
716  U[2][6] = (1-K[4])*K[5];
717  U[2][7] = K[4]*K[5];
718  }
719  else if (nvert == 4)
720  {
721  // linear tetr
722  // for regular reference tetrahedron
723  U[0][0] = -1;
724  U[0][1] = 1;
725  U[0][2] = 0;
726  U[0][3] = 0;
727 
728  U[1][0] = -1./sqrt3;
729  U[1][1] = -1./sqrt3;
730  U[1][2] = 2./sqrt3;
731  U[1][3] = 0;
732 
733  U[2][0] = -1./sqrt6;
734  U[2][1] = -1./sqrt6;
735  U[2][2] = -1./sqrt6;
736  U[2][3] = 3./sqrt6;
737 
738  // for right corner reference tetrahedron
739  // U[0][0] = -1; U[1][0] = -1; U[2][0] = -1;
740  // U[0][1] = 1; U[1][1] = 0; U[2][1] = 0;
741  // U[0][2] = 0; U[1][2] = 1; U[2][2] = 0;
742  // U[0][3] = 0; U[1][3] = 0; U[2][3] = 1;
743  }
744  else if (nvert == 6)
745  {
746  // prism
747  // for regular triangle in the prism base
748  U[0][0] = -1+K[0];
749  U[0][1] = 1-K[0];
750  U[0][2] = 0;
751  U[0][3] = -K[0];
752  U[0][4] = K[0];
753  U[0][5] = 0;
754 
755  U[1][0] = 0.5*(-1+K[1]);
756  U[1][1] = 0.5*(-1+K[1]);
757  U[1][2] = (1-K[1]);
758  U[1][3] = -0.5*K[1];
759  U[1][4] = -0.5*K[1];
760  U[1][5] = K[1];
761 
762  U[2][0] = -1. + K[2] + 0.5*K[3];
763  U[2][1] = -K[2] + 0.5*K[3];
764  U[2][2] = -K[3];
765  U[2][3] = 1 - K[2] - 0.5*K[3];
766  U[2][4] = K[2] - 0.5*K[3];
767  U[2][5] = K[3];
768  }
769  else if (nvert == 10)
770  {
771  // quad tetr
772  U[0][0] = -3.*(1 - K[0] - 0.5*K[1] - K[2]/3.) + (K[0] - 0.5*K[1] - K[2]/3.) + (K[1] - K[2]/3.) + K[2];
773  U[0][1] = -1.*(1 - K[0] - 0.5*K[1] - K[2]/3.) + 3.*(K[0] - 0.5*K[1] - K[2]/3.) - (K[1] - K[2]/3.) - K[2];
774  U[0][2] = 0.;
775  U[0][3] = 0.;
776  U[0][4] = 4.*(K[1] - K[2]/3.);
777  U[0][5] = -4.*(K[1] - K[2]/3.);
778  U[0][6] = 4.*(1. - K[0] - 0.5*K[1] - K[2]/3.) - 4.*(K[0] - 0.5*K[1] - K[2]/3.);
779  U[0][7] = 4.*K[2];
780  U[0][8] = 0.;
781  U[0][9] = -4.*K[2];
782 
783  U[1][0] = -1.5*(1. - K[3] - K[4]*0.5 - K[5]/3.) + 0.5*(K[3] - 0.5*K[4] - K[5]/3.) + 0.5*(K[4] - K[5]/3.) + 0.5*K[5];
784  U[1][1] = 0.5*(1. - K[3] - K[4]*0.5 - K[5]/3.) - 1.5*(K[3] - 0.5*K[4] - K[5]/3.) + 0.5*(K[4] - K[5]/3.) + 0.5*K[5];
785  U[1][2] = -1.*(1. - K[3] - K[4]*0.5 - K[5]/3.) - (K[3] - 0.5*K[4] - K[5]/3.) + 3.*(K[4] - K[5]/3.) - K[5];
786  U[1][3] = 0.;
787  U[1][4] = 4.*( K[3] - 0.5*K[4] - K[5]/3.) - 2.*(K[4] - K[5]/3.);
788  U[1][5] = 4.*(1. - K[3] - 0.5*K[4] - K[5]/3.) - 2.*(K[4] - K[5]/3.);
789  U[1][6] = -2.*(1. - K[3] - 0.5*K[4] - K[5]/3.) - 2.*(K[3] - 0.5*K[4] - K[5]/3.);
790  U[1][7] = -2.*K[5];
791  U[1][8] = 4.*K[5];
792  U[1][9] = -2.*K[5];
793 
794  U[2][0] = -(1. - K[6] - 0.5*K[7] - K[8]/3.) + (K[6] - 0.5*K[7] - K[8]/3.)/3. + (K[7] - K[8]/3.)/3. + K[8]/3.;
795  U[2][1] = (1. - K[6] - 0.5*K[7] - K[8]/3.)/3. - (K[6] - 0.5*K[7] - K[8]/3.) + (K[7] - K[8]/3.)/3. + K[8]/3.;
796  U[2][2] = (1. - K[6] - 0.5*K[7] - K[8]/3.)/3. + (K[6] - 0.5*K[7] - K[8]/3.)/3. - (K[7] - K[8]/3.) + K[8]/3.;
797  U[2][3] = -(1. - K[6] - 0.5*K[7] - K[8]/3.) - (K[6] - 0.5*K[7] - K[8]/3.) - (K[7] - K[8]/3.) + 3.*K[8];
798  U[2][4] = -4.*(K[6] - K[7]*0.5 - K[8]/3.)/3. - 4.*(K[7] - K[8]/3.)/3.;
799  U[2][5] = -4.*(1. - K[6] - K[7]*0.5 - K[8]/3.)/3. - 4.*( K[7] - K[8]/3.)/3.;
800  U[2][6] = -4.*(1. - K[6] - K[7]*0.5 - K[8]/3.)/3. - 4.*(K[6] - 0.5*K[7] - K[8]/3.)/3.;
801  U[2][7] = 4.*(K[6] - K[7]*0.5 - K[8]/3.) - 4.*K[8]/3.;
802  U[2][8] = 4.*(K[7] - K[8]/3.) - 4.*K[8]/3.;
803  U[2][9] = 4.*(1. - K[6] - K[7]*0.5 - K[8]/3.) - 4.*K[8]/3.;
804  }
805  else
806  libmesh_error_msg("Invalid nvert = " << nvert);
807  }
808 
809  if (me == 1)
810  Q = U;
811 
812  else
813  {
814  for (unsigned i=0; i<_dim; i++)
815  for (int j=0; j<nvert; j++)
816  {
817  Q[i][j] = 0;
818  for (unsigned k=0; k<_dim; k++)
819  Q[i][j] += H[k][i]*U[k][j];
820  }
821  }
822 
823  return 0;
824 }
825 
826 
827 
828 // Specify adaptive function
829 void VariationalMeshSmoother::adp_renew(const Array2D<double> & R,
830  const Array2D<int> & cells,
831  std::vector<double> & afun,
832  int adp)
833 {
834  // evaluates given adaptive function on the provided mesh
835  if (adp < 0)
836  {
837  for (dof_id_type i=0; i<_n_cells; i++)
838  {
839  double
840  x = 0.,
841  y = 0.,
842  z = 0.;
843  int nvert = 0;
844 
845  while (cells[i][nvert] >= 0)
846  {
847  x += R[cells[i][nvert]][0];
848  y += R[cells[i][nvert]][1];
849  if (_dim == 3)
850  z += R[cells[i][nvert]][2];
851  nvert++;
852  }
853 
854  // adaptive function, cell based
855  afun[i] = 5.*(x+y+z);
856  }
857  }
858 
859  else if (adp > 0)
860  {
861  for (dof_id_type i=0; i<_n_nodes; i++)
862  {
863  double z = 0;
864  for (unsigned j=0; j<_dim; j++)
865  z += R[i][j];
866 
867  // adaptive function, node based
868  afun[i] = 5*std::sin(R[i][0]);
869  }
870  }
871 }
872 
873 
874 
875 // Preprocess mesh data and control smoothing/untangling iterations
876 void VariationalMeshSmoother::full_smooth(Array2D<double> & R,
877  const std::vector<int> & mask,
878  const Array2D<int> & cells,
879  const std::vector<int> & mcells,
880  const std::vector<int> & edges,
881  const std::vector<int> & hnodes,
882  double w,
883  const std::vector<int> & iter,
884  int me,
885  const Array3D<double> & H,
886  int adp,
887  int gr)
888 {
889  // Control the amount of print statements in this function
890  int msglev = 1;
891 
892  dof_id_type afun_size = 0;
893 
894  // Adaptive function is on cells
895  if (adp < 0)
896  afun_size = _n_cells;
897 
898  // Adaptive function is on nodes
899  else if (adp > 0)
900  afun_size = _n_nodes;
901 
902  std::vector<double> afun(afun_size);
903  std::vector<int> maskf(_n_nodes);
904  std::vector<double> Gamma(_n_cells);
905 
906  if (msglev >= 1)
907  _logfile << "N=" << _n_nodes
908  << " ncells=" << _n_cells
909  << " nedges=" << _n_hanging_edges
910  << std::endl;
911 
912 
913  // Boundary node counting
914  int NBN=0;
915  for (dof_id_type i=0; i<_n_nodes; i++)
916  if (mask[i] == 2 || mask[i] == 1)
917  NBN++;
918 
919  if (NBN > 0)
920  {
921  if (msglev >= 1)
922  _logfile << "# of Boundary Nodes=" << NBN << std::endl;
923 
924  NBN=0;
925  for (dof_id_type i=0; i<_n_nodes; i++)
926  if (mask[i] == 2)
927  NBN++;
928 
929  if (msglev >= 1)
930  _logfile << "# of moving Boundary Nodes=" << NBN << std::endl;
931  }
932 
933  for (dof_id_type i=0; i<_n_nodes; i++)
934  {
935  if ((mask[i]==1) || (mask[i]==2))
936  maskf[i] = 1;
937  else
938  maskf[i] = 0;
939  }
940 
941  // determination of min jacobian
942  double vol, Vmin;
943  double qmin = minq(R, cells, mcells, me, H, vol, Vmin);
944 
945  if (me > 1)
946  vol = 1.;
947 
948  if (msglev >= 1)
949  _logfile << "vol=" << vol
950  << " qmin=" << qmin
951  << " min volume = " << Vmin
952  << std::endl;
953 
954  // compute max distortion measure over all cells
955  double epsilon = 1.e-9;
956  double eps = qmin < 0 ? std::sqrt(epsilon*epsilon+0.004*qmin*qmin*vol*vol) : epsilon;
957  double emax = maxE(R, cells, mcells, me, H, vol, eps, w, Gamma, qmin);
958 
959  if (msglev >= 1)
960  _logfile << " emax=" << emax << std::endl;
961 
962  // unfolding/smoothing
963 
964  // iteration tolerance
965  double Enm1 = 1.;
966 
967  // read adaptive function from file
968  if (adp*gr != 0)
969  read_adp(afun);
970 
971  {
972  int
973  counter = 0,
974  ii = 0;
975 
976  while (((qmin <= 0) || (counter < iter[0]) || (std::abs(emax-Enm1) > 1e-3)) &&
977  (ii < iter[1]) &&
978  (counter < iter[1]))
979  {
980  libmesh_assert_less (counter, iter[1]);
981 
982  Enm1 = emax;
983 
984  if ((ii >= 0) && (ii % 2 == 0))
985  {
986  if (qmin < 0)
987  eps = std::sqrt(epsilon*epsilon + 0.004*qmin*qmin*vol*vol);
988  else
989  eps = epsilon;
990  }
991 
992  int ladp = adp;
993 
994  if ((qmin <= 0) || (counter < ii))
995  ladp = 0;
996 
997  // update adaptation function before each iteration
998  if ((ladp != 0) && (gr == 0))
999  adp_renew(R, cells, afun, adp);
1000 
1001  double Jk = minJ(R, maskf, cells, mcells, eps, w, me, H, vol, edges, hnodes,
1002  msglev, Vmin, emax, qmin, ladp, afun);
1003 
1004  if (qmin > 0)
1005  counter++;
1006  else
1007  ii++;
1008 
1009  if (msglev >= 1)
1010  _logfile << "niter=" << counter
1011  << ", qmin*G/vol=" << qmin
1012  << ", Vmin=" << Vmin
1013  << ", emax=" << emax
1014  << ", Jk=" << Jk
1015  << ", Enm1=" << Enm1
1016  << std::endl;
1017  }
1018  }
1019 
1020  // BN correction - 2D only!
1021  epsilon = 1.e-9;
1022  if (NBN > 0)
1023  for (int counter=0; counter<iter[2]; counter++)
1024  {
1025  // update adaptation function before each iteration
1026  if ((adp != 0) && (gr == 0))
1027  adp_renew(R, cells, afun, adp);
1028 
1029  double Jk = minJ_BC(R, mask, cells, mcells, eps, w, me, H, vol, msglev, Vmin, emax, qmin, adp, afun, NBN);
1030 
1031  if (msglev >= 1)
1032  _logfile << "NBC niter=" << counter
1033  << ", qmin*G/vol=" << qmin
1034  << ", Vmin=" << Vmin
1035  << ", emax=" << emax
1036  << std::endl;
1037 
1038  // Outrageous Enm1 to make sure we hit this at least once
1039  Enm1 = 99999;
1040 
1041  // Now that we've moved the boundary nodes (or not) we need to resmooth
1042  for (int j=0; j<iter[1]; j++)
1043  {
1044  if (std::abs(emax-Enm1) < 1e-2)
1045  break;
1046 
1047  // Save off the error from the previous smoothing step
1048  Enm1 = emax;
1049 
1050  // update adaptation function before each iteration
1051  if ((adp != 0) && (gr == 0))
1052  adp_renew(R, cells, afun, adp);
1053 
1054  Jk = minJ(R, maskf, cells, mcells, eps, w, me, H, vol, edges, hnodes, msglev, Vmin, emax, qmin, adp, afun);
1055 
1056  if (msglev >= 1)
1057  _logfile << " Re-smooth: niter=" << j
1058  << ", qmin*G/vol=" << qmin
1059  << ", Vmin=" << Vmin
1060  << ", emax=" << emax
1061  << ", Jk=" << Jk
1062  << std::endl;
1063  }
1064 
1065  if (msglev >= 1)
1066  _logfile << "NBC smoothed niter=" << counter
1067  << ", qmin*G/vol=" << qmin
1068  << ", Vmin=" << Vmin
1069  << ", emax=" << emax
1070  << std::endl;
1071  }
1072 }
1073 
1074 
1075 
1076 // Determines the values of maxE_theta
1077 double VariationalMeshSmoother::maxE(Array2D<double> & R,
1078  const Array2D<int> & cells,
1079  const std::vector<int> & mcells,
1080  int me,
1081  const Array3D<double> & H,
1082  double v,
1083  double epsilon,
1084  double w,
1085  std::vector<double> & Gamma,
1086  double & qmin)
1087 {
1088  Array2D<double> Q(3, 3*_dim + _dim%2);
1089  std::vector<double> K(9);
1090 
1091  double
1092  gemax = -1.e32,
1093  vmin = 1.e32;
1094 
1095  for (dof_id_type ii=0; ii<_n_cells; ii++)
1096  if (mcells[ii] >= 0)
1097  {
1098  // Final value of E will be saved in Gamma at the end of this loop
1099  double E = 0.;
1100 
1101  if (_dim == 2)
1102  {
1103  if (cells[ii][3] == -1)
1104  {
1105  // tri
1106  basisA(Q, 3, K, H[ii], me);
1107 
1108  std::vector<double> a1(3), a2(3);
1109  for (int k=0; k<2; k++)
1110  for (int l=0; l<3; l++)
1111  {
1112  a1[k] += Q[k][l]*R[cells[ii][l]][0];
1113  a2[k] += Q[k][l]*R[cells[ii][l]][1];
1114  }
1115 
1116  double det = jac2(a1[0], a1[1], a2[0], a2[1]);
1117  double tr = 0.5*(a1[0]*a1[0] + a2[0]*a2[0] + a1[1]*a1[1] + a2[1]*a2[1]);
1118  double chi = 0.5*(det+std::sqrt(det*det+epsilon*epsilon));
1119  E = (1-w)*tr/chi + 0.5*w*(v + det*det/v)/chi;
1120 
1121  if (E > gemax)
1122  gemax = E;
1123  if (vmin > det)
1124  vmin = det;
1125 
1126  }
1127  else if (cells[ii][4] == -1)
1128  {
1129  // quad
1130  for (int i=0; i<2; i++)
1131  {
1132  K[0] = i;
1133  for (int j=0; j<2; j++)
1134  {
1135  K[1] = j;
1136  basisA(Q, 4, K, H[ii], me);
1137 
1138  std::vector<double> a1(3), a2(3);
1139  for (int k=0; k<2; k++)
1140  for (int l=0; l<4; l++)
1141  {
1142  a1[k] += Q[k][l]*R[cells[ii][l]][0];
1143  a2[k] += Q[k][l]*R[cells[ii][l]][1];
1144  }
1145 
1146  double det = jac2(a1[0],a1[1],a2[0],a2[1]);
1147  double tr = 0.5*(a1[0]*a1[0] + a2[0]*a2[0] + a1[1]*a1[1] + a2[1]*a2[1]);
1148  double chi = 0.5*(det+std::sqrt(det*det+epsilon*epsilon));
1149  E += 0.25*((1-w)*tr/chi + 0.5*w*(v + det*det/v)/chi);
1150 
1151  if (vmin > det)
1152  vmin = det;
1153  }
1154  }
1155 
1156  if (E > gemax)
1157  gemax = E;
1158  }
1159  else
1160  {
1161  // quad tri
1162  for (int i=0; i<3; i++)
1163  {
1164  K[0] = i*0.5;
1165  int k = i/2;
1166  K[1] = static_cast<double>(k);
1167 
1168  for (int j=0; j<3; j++)
1169  {
1170  K[2] = j*0.5;
1171  k = j/2;
1172  K[3] = static_cast<double>(k);
1173 
1174  basisA(Q, 6, K, H[ii], me);
1175 
1176  std::vector<double> a1(3), a2(3);
1177  for (int k2=0; k2<2; k2++)
1178  for (int l=0; l<6; l++)
1179  {
1180  a1[k2] += Q[k2][l]*R[cells[ii][l]][0];
1181  a2[k2] += Q[k2][l]*R[cells[ii][l]][1];
1182  }
1183 
1184  double det = jac2(a1[0],a1[1],a2[0],a2[1]);
1185  double sigma = 1./24.;
1186 
1187  if (i==j)
1188  sigma = 1./12.;
1189 
1190  double tr = 0.5*(a1[0]*a1[0] + a2[0]*a2[0] + a1[1]*a1[1] + a2[1]*a2[1]);
1191  double chi = 0.5*(det + std::sqrt(det*det + epsilon*epsilon));
1192  E += sigma*((1-w)*tr/chi + 0.5*w*(v + det*det/v)/chi);
1193  if (vmin > det)
1194  vmin = det;
1195  }
1196  }
1197 
1198  if (E > gemax)
1199  gemax = E;
1200  }
1201  }
1202 
1203  if (_dim == 3)
1204  {
1205  if (cells[ii][4] == -1)
1206  {
1207  // tetr
1208  basisA(Q, 4, K, H[ii], me);
1209 
1210  std::vector<double> a1(3), a2(3), a3(3);
1211  for (int k=0; k<3; k++)
1212  for (int l=0; l<4; l++)
1213  {
1214  a1[k] += Q[k][l]*R[cells[ii][l]][0];
1215  a2[k] += Q[k][l]*R[cells[ii][l]][1];
1216  a3[k] += Q[k][l]*R[cells[ii][l]][2];
1217  }
1218 
1219  double det = jac3(a1[0], a1[1], a1[2],
1220  a2[0], a2[1], a2[2],
1221  a3[0], a3[1], a3[2]);
1222  double tr = 0.;
1223  for (int k=0; k<3; k++)
1224  tr += (a1[k]*a1[k] + a2[k]*a2[k] + a3[k]*a3[k])/3.;
1225 
1226  double chi = 0.5*(det+std::sqrt(det*det+epsilon*epsilon));
1227  E = (1-w)*pow(tr,1.5)/chi + 0.5*w*(v + det*det/v)/chi;
1228 
1229  if (E > gemax)
1230  gemax = E;
1231 
1232  if (vmin > det)
1233  vmin = det;
1234  }
1235  else if (cells[ii][6] == -1)
1236  {
1237  // prism
1238  for (int i=0; i<2; i++)
1239  {
1240  K[0] = i;
1241  for (int j=0; j<2; j++)
1242  {
1243  K[1] = j;
1244  for (int k=0; k<3; k++)
1245  {
1246  K[2] = 0.5*static_cast<double>(k);
1247  K[3] = static_cast<double>(k % 2);
1248  basisA(Q, 6, K, H[ii], me);
1249 
1250  std::vector<double> a1(3), a2(3), a3(3);
1251  for (int kk=0; kk<3; kk++)
1252  for (int ll=0; ll<6; ll++)
1253  {
1254  a1[kk] += Q[kk][ll]*R[cells[ii][ll]][0];
1255  a2[kk] += Q[kk][ll]*R[cells[ii][ll]][1];
1256  a3[kk] += Q[kk][ll]*R[cells[ii][ll]][2];
1257  }
1258 
1259  double det = jac3(a1[0], a1[1], a1[2],
1260  a2[0], a2[1], a2[2],
1261  a3[0], a3[1], a3[2]);
1262  double tr = 0;
1263  for (int kk=0; kk<3; kk++)
1264  tr += (a1[kk]*a1[kk] + a2[kk]*a2[kk] + a3[kk]*a3[kk])/3.;
1265 
1266  double chi = 0.5*(det+std::sqrt(det*det+epsilon*epsilon));
1267  E += ((1-w)*pow(tr,1.5)/chi + 0.5*w*(v + det*det/v)/chi)/12.;
1268  if (vmin > det)
1269  vmin = det;
1270  }
1271  }
1272  }
1273 
1274  if (E > gemax)
1275  gemax = E;
1276  }
1277  else if (cells[ii][8] == -1)
1278  {
1279  // hex
1280  for (int i=0; i<2; i++)
1281  {
1282  K[0] = i;
1283  for (int j=0; j<2; j++)
1284  {
1285  K[1] = j;
1286  for (int k=0; k<2; k++)
1287  {
1288  K[2] = k;
1289  for (int l=0; l<2; l++)
1290  {
1291  K[3] = l;
1292  for (int m=0; m<2; m++)
1293  {
1294  K[4] = m;
1295  for (int nn=0; nn<2; nn++)
1296  {
1297  K[5] = nn;
1298  basisA(Q, 8, K, H[ii], me);
1299 
1300  std::vector<double> a1(3), a2(3), a3(3);
1301  for (int kk=0; kk<3; kk++)
1302  for (int ll=0; ll<8; ll++)
1303  {
1304  a1[kk] += Q[kk][ll]*R[cells[ii][ll]][0];
1305  a2[kk] += Q[kk][ll]*R[cells[ii][ll]][1];
1306  a3[kk] += Q[kk][ll]*R[cells[ii][ll]][2];
1307  }
1308 
1309  double det = jac3(a1[0], a1[1], a1[2],
1310  a2[0], a2[1], a2[2],
1311  a3[0], a3[1], a3[2]);
1312  double sigma = 0.;
1313 
1314  if ((i==nn) && (j==l) && (k==m))
1315  sigma = 1./27.;
1316 
1317  if (((i==nn) && (j==l) && (k!=m)) ||
1318  ((i==nn) && (j!=l) && (k==m)) ||
1319  ((i!=nn) && (j==l) && (k==m)))
1320  sigma = 1./54.;
1321 
1322  if (((i==nn) && (j!=l) && (k!=m)) ||
1323  ((i!=nn) && (j!=l) && (k==m)) ||
1324  ((i!=nn) && (j==l) && (k!=m)))
1325  sigma = 1./108.;
1326 
1327  if ((i!=nn) && (j!=l) && (k!=m))
1328  sigma = 1./216.;
1329 
1330  double tr = 0;
1331  for (int kk=0; kk<3; kk++)
1332  tr += (a1[kk]*a1[kk] + a2[kk]*a2[kk] + a3[kk]*a3[kk])/3.;
1333 
1334  double chi = 0.5*(det+std::sqrt(det*det + epsilon*epsilon));
1335  E += ((1-w)*pow(tr,1.5)/chi + 0.5*w*(v + det*det/v)/chi)*sigma;
1336  if (vmin > det)
1337  vmin = det;
1338  }
1339  }
1340  }
1341  }
1342  }
1343  }
1344  if (E > gemax)
1345  gemax = E;
1346  }
1347  else
1348  {
1349  // quad tetr
1350  for (int i=0; i<4; i++)
1351  {
1352  for (int j=0; j<4; j++)
1353  {
1354  for (int k=0; k<4; k++)
1355  {
1356  switch (i)
1357  {
1358  case 0:
1359  K[0] = 0;
1360  K[1] = 0;
1361  K[2] = 0;
1362  break;
1363  case 1:
1364  K[0] = 1;
1365  K[1] = 0;
1366  K[2] = 0;
1367  break;
1368  case 2:
1369  K[0] = 0.5;
1370  K[1] = 1;
1371  K[2] = 0;
1372  break;
1373  case 3:
1374  K[0] = 0.5;
1375  K[1] = 1./3.;
1376  K[2] = 1;
1377  break;
1378  default:
1379  break;
1380  }
1381 
1382  switch (j)
1383  {
1384  case 0:
1385  K[3] = 0;
1386  K[4] = 0;
1387  K[5] = 0;
1388  break;
1389  case 1:
1390  K[3] = 1;
1391  K[4] = 0;
1392  K[5] = 0;
1393  break;
1394  case 2:
1395  K[3] = 0.5;
1396  K[4] = 1;
1397  K[5] = 0;
1398  break;
1399  case 3:
1400  K[3] = 0.5;
1401  K[4] = 1./3.;
1402  K[5] = 1;
1403  break;
1404  default:
1405  break;
1406  }
1407 
1408  switch (k)
1409  {
1410  case 0:
1411  K[6] = 0;
1412  K[7] = 0;
1413  K[8] = 0;
1414  break;
1415  case 1:
1416  K[6] = 1;
1417  K[7] = 0;
1418  K[8] = 0;
1419  break;
1420  case 2:
1421  K[6] = 0.5;
1422  K[7] = 1;
1423  K[8] = 0;
1424  break;
1425  case 3:
1426  K[6] = 0.5;
1427  K[7] = 1./3.;
1428  K[8] = 1;
1429  break;
1430  default:
1431  break;
1432  }
1433 
1434  basisA(Q, 10, K, H[ii], me);
1435 
1436  std::vector<double> a1(3), a2(3), a3(3);
1437  for (int kk=0; kk<3; kk++)
1438  for (int ll=0; ll<10; ll++)
1439  {
1440  a1[kk] += Q[kk][ll]*R[cells[ii][ll]][0];
1441  a2[kk] += Q[kk][ll]*R[cells[ii][ll]][1];
1442  a3[kk] += Q[kk][ll]*R[cells[ii][ll]][2];
1443  }
1444 
1445  double det = jac3(a1[0], a1[1], a1[2],
1446  a2[0], a2[1], a2[2],
1447  a3[0], a3[1], a3[2]);
1448  double sigma = 0.;
1449 
1450  if ((i==j) && (j==k))
1451  sigma = 1./120.;
1452  else if ((i==j) || (j==k) || (i==k))
1453  sigma = 1./360.;
1454  else
1455  sigma = 1./720.;
1456 
1457  double tr = 0;
1458  for (int kk=0; kk<3; kk++)
1459  tr += (a1[kk]*a1[kk] + a2[kk]*a2[kk] + a3[kk]*a3[kk])/3.;
1460 
1461  double chi = 0.5*(det+std::sqrt(det*det+epsilon*epsilon));
1462  E += ((1-w)*pow(tr,1.5)/chi + 0.5*w*(v+det*det/v)/chi)*sigma;
1463  if (vmin > det)
1464  vmin = det;
1465  }
1466  }
1467  }
1468 
1469  if (E > gemax)
1470  gemax = E;
1471  }
1472  }
1473  Gamma[ii] = E;
1474  }
1475 
1476  qmin = vmin;
1477 
1478  return gemax;
1479 }
1480 
1481 
1482 
1483 // Compute min Jacobian determinant (minq), min cell volume (Vmin), and average cell volume (vol).
1484 double VariationalMeshSmoother::minq(const Array2D<double> & R,
1485  const Array2D<int> & cells,
1486  const std::vector<int> & mcells,
1487  int me,
1488  const Array3D<double> & H,
1489  double & vol,
1490  double & Vmin)
1491 {
1492  std::vector<double> K(9);
1493  Array2D<double> Q(3, 3*_dim + _dim%2);
1494 
1495  double v = 0;
1496  double vmin = 1.e32;
1497  double gqmin = 1.e32;
1498 
1499  for (dof_id_type ii=0; ii<_n_cells; ii++)
1500  if (mcells[ii] >= 0)
1501  {
1502  if (_dim == 2)
1503  {
1504  // 2D
1505  if (cells[ii][3] == -1)
1506  {
1507  // tri
1508  basisA(Q, 3, K, H[ii], me);
1509 
1510  std::vector<double> a1(3), a2(3);
1511  for (int k=0; k<2; k++)
1512  for (int l=0; l<3; l++)
1513  {
1514  a1[k] += Q[k][l]*R[cells[ii][l]][0];
1515  a2[k] += Q[k][l]*R[cells[ii][l]][1];
1516  }
1517 
1518  double det = jac2(a1[0],a1[1],a2[0],a2[1]);
1519  if (gqmin > det)
1520  gqmin = det;
1521 
1522  if (vmin > det)
1523  vmin = det;
1524 
1525  v += det;
1526  }
1527  else if (cells[ii][4] == -1)
1528  {
1529  // quad
1530  double vcell = 0.;
1531  for (int i=0; i<2; i++)
1532  {
1533  K[0] = i;
1534  for (int j=0; j<2; j++)
1535  {
1536  K[1] = j;
1537  basisA(Q, 4, K, H[ii], me);
1538 
1539  std::vector<double> a1(3), a2(3);
1540  for (int k=0; k<2; k++)
1541  for (int l=0; l<4; l++)
1542  {
1543  a1[k] += Q[k][l]*R[cells[ii][l]][0];
1544  a2[k] += Q[k][l]*R[cells[ii][l]][1];
1545  }
1546 
1547  double det = jac2(a1[0],a1[1],a2[0],a2[1]);
1548  if (gqmin > det)
1549  gqmin = det;
1550 
1551  v += 0.25*det;
1552  vcell += 0.25*det;
1553  }
1554  }
1555  if (vmin > vcell)
1556  vmin = vcell;
1557  }
1558  else
1559  {
1560  // quad tri
1561  double vcell = 0.;
1562  for (int i=0; i<3; i++)
1563  {
1564  K[0] = i*0.5;
1565  int k = i/2;
1566  K[1] = static_cast<double>(k);
1567 
1568  for (int j=0; j<3; j++)
1569  {
1570  K[2] = j*0.5;
1571  k = j/2;
1572  K[3] = static_cast<double>(k);
1573  basisA(Q, 6, K, H[ii], me);
1574 
1575  std::vector<double> a1(3), a2(3);
1576  for (int k2=0; k2<2; k2++)
1577  for (int l=0; l<6; l++)
1578  {
1579  a1[k2] += Q[k2][l]*R[cells[ii][l]][0];
1580  a2[k2] += Q[k2][l]*R[cells[ii][l]][1];
1581  }
1582 
1583  double det = jac2(a1[0], a1[1], a2[0], a2[1]);
1584  if (gqmin > det)
1585  gqmin = det;
1586 
1587  double sigma = 1./24.;
1588  if (i == j)
1589  sigma = 1./12.;
1590 
1591  v += sigma*det;
1592  vcell += sigma*det;
1593  }
1594  }
1595  if (vmin > vcell)
1596  vmin = vcell;
1597  }
1598  }
1599  if (_dim == 3)
1600  {
1601  // 3D
1602  if (cells[ii][4] == -1)
1603  {
1604  // tetr
1605  basisA(Q, 4, K, H[ii], me);
1606 
1607  std::vector<double> a1(3), a2(3), a3(3);
1608  for (int k=0; k<3; k++)
1609  for (int l=0; l<4; l++)
1610  {
1611  a1[k] += Q[k][l]*R[cells[ii][l]][0];
1612  a2[k] += Q[k][l]*R[cells[ii][l]][1];
1613  a3[k] += Q[k][l]*R[cells[ii][l]][2];
1614  }
1615 
1616  double det = jac3(a1[0], a1[1], a1[2],
1617  a2[0], a2[1], a2[2],
1618  a3[0], a3[1], a3[2]);
1619 
1620  if (gqmin > det)
1621  gqmin = det;
1622 
1623  if (vmin > det)
1624  vmin = det;
1625  v += det;
1626  }
1627  else if (cells[ii][6] == -1)
1628  {
1629  // prism
1630  double vcell = 0.;
1631  for (int i=0; i<2; i++)
1632  {
1633  K[0] = i;
1634  for (int j=0; j<2; j++)
1635  {
1636  K[1] = j;
1637 
1638  for (int k=0; k<3; k++)
1639  {
1640  K[2] = 0.5*k;
1641  K[3] = static_cast<double>(k%2);
1642  basisA(Q, 6, K, H[ii], me);
1643 
1644  std::vector<double> a1(3), a2(3), a3(3);
1645  for (int kk=0; kk<3; kk++)
1646  for (int ll=0; ll<6; ll++)
1647  {
1648  a1[kk] += Q[kk][ll]*R[cells[ii][ll]][0];
1649  a2[kk] += Q[kk][ll]*R[cells[ii][ll]][1];
1650  a3[kk] += Q[kk][ll]*R[cells[ii][ll]][2];
1651  }
1652 
1653  double det = jac3(a1[0], a1[1], a1[2],
1654  a2[0], a2[1], a2[2],
1655  a3[0], a3[1], a3[2]);
1656  if (gqmin > det)
1657  gqmin = det;
1658 
1659  double sigma = 1./12.;
1660  v += sigma*det;
1661  vcell += sigma*det;
1662  }
1663  }
1664  }
1665  if (vmin > vcell)
1666  vmin = vcell;
1667  }
1668  else if (cells[ii][8] == -1)
1669  {
1670  // hex
1671  double vcell = 0.;
1672  for (int i=0; i<2; i++)
1673  {
1674  K[0] = i;
1675  for (int j=0; j<2; j++)
1676  {
1677  K[1] = j;
1678  for (int k=0; k<2; k++)
1679  {
1680  K[2] = k;
1681  for (int l=0; l<2; l++)
1682  {
1683  K[3] = l;
1684  for (int m=0; m<2; m++)
1685  {
1686  K[4] = m;
1687  for (int nn=0; nn<2; nn++)
1688  {
1689  K[5] = nn;
1690  basisA(Q, 8, K, H[ii], me);
1691 
1692  std::vector<double> a1(3), a2(3), a3(3);
1693  for (int kk=0; kk<3; kk++)
1694  for (int ll=0; ll<8; ll++)
1695  {
1696  a1[kk] += Q[kk][ll]*R[cells[ii][ll]][0];
1697  a2[kk] += Q[kk][ll]*R[cells[ii][ll]][1];
1698  a3[kk] += Q[kk][ll]*R[cells[ii][ll]][2];
1699  }
1700 
1701  double det = jac3(a1[0], a1[1], a1[2],
1702  a2[0], a2[1], a2[2],
1703  a3[0], a3[1], a3[2]);
1704 
1705  if (gqmin > det)
1706  gqmin = det;
1707 
1708  double sigma = 0.;
1709 
1710  if ((i==nn) && (j==l) && (k==m))
1711  sigma = 1./27.;
1712 
1713  if (((i==nn) && (j==l) && (k!=m)) ||
1714  ((i==nn) && (j!=l) && (k==m)) ||
1715  ((i!=nn) && (j==l) && (k==m)))
1716  sigma = 1./54.;
1717 
1718  if (((i==nn) && (j!=l) && (k!=m)) ||
1719  ((i!=nn) && (j!=l) && (k==m)) ||
1720  ((i!=nn) && (j==l) && (k!=m)))
1721  sigma = 1./108.;
1722 
1723  if ((i!=nn) && (j!=l) && (k!=m))
1724  sigma = 1./216.;
1725 
1726  v += sigma*det;
1727  vcell += sigma*det;
1728  }
1729  }
1730  }
1731  }
1732  }
1733  }
1734 
1735  if (vmin > vcell)
1736  vmin = vcell;
1737  }
1738  else
1739  {
1740  // quad tetr
1741  double vcell = 0.;
1742  for (int i=0; i<4; i++)
1743  {
1744  for (int j=0; j<4; j++)
1745  {
1746  for (int k=0; k<4; k++)
1747  {
1748  switch (i)
1749  {
1750  case 0:
1751  K[0] = 0;
1752  K[1] = 0;
1753  K[2] = 0;
1754  break;
1755 
1756  case 1:
1757  K[0] = 1;
1758  K[1] = 0;
1759  K[2] = 0;
1760  break;
1761 
1762  case 2:
1763  K[0] = 0.5;
1764  K[1] = 1;
1765  K[2] = 0;
1766  break;
1767 
1768  case 3:
1769  K[0] = 0.5;
1770  K[1] = 1./3.;
1771  K[2] = 1;
1772  break;
1773 
1774  default:
1775  break;
1776  }
1777  switch (j)
1778  {
1779  case 0:
1780  K[3] = 0;
1781  K[4] = 0;
1782  K[5] = 0;
1783  break;
1784 
1785  case 1:
1786  K[3] = 1;
1787  K[4] = 0;
1788  K[5] = 0;
1789  break;
1790 
1791  case 2:
1792  K[3] = 0.5;
1793  K[4] = 1;
1794  K[5] = 0;
1795  break;
1796 
1797  case 3:
1798  K[3] = 0.5;
1799  K[4] = 1./3.;
1800  K[5] = 1;
1801  break;
1802 
1803  default:
1804  break;
1805  }
1806  switch (k)
1807  {
1808  case 0:
1809  K[6] = 0;
1810  K[7] = 0;
1811  K[8] = 0;
1812  break;
1813 
1814  case 1:
1815  K[6] = 1;
1816  K[7] = 0;
1817  K[8] = 0;
1818  break;
1819 
1820  case 2:
1821  K[6] = 0.5;
1822  K[7] = 1;
1823  K[8] = 0;
1824  break;
1825 
1826  case 3:
1827  K[6] = 0.5;
1828  K[7] = 1./3.;
1829  K[8] = 1;
1830  break;
1831 
1832  default:
1833  break;
1834  }
1835 
1836  basisA(Q, 10, K, H[ii], me);
1837 
1838  std::vector<double> a1(3), a2(3), a3(3);
1839  for (int kk=0; kk<3; kk++)
1840  for (int ll=0; ll<10; ll++)
1841  {
1842  a1[kk] += Q[kk][ll]*R[cells[ii][ll]][0];
1843  a2[kk] += Q[kk][ll]*R[cells[ii][ll]][1];
1844  a3[kk] += Q[kk][ll]*R[cells[ii][ll]][2];
1845  }
1846 
1847  double det = jac3(a1[0], a1[1], a1[2],
1848  a2[0], a2[1], a2[2],
1849  a3[0], a3[1], a3[2]);
1850  if (gqmin > det)
1851  gqmin = det;
1852 
1853  double sigma = 0.;
1854 
1855  if ((i==j) && (j==k))
1856  sigma = 1./120.;
1857 
1858  else if ((i==j) || (j==k) || (i==k))
1859  sigma = 1./360.;
1860 
1861  else
1862  sigma = 1./720.;
1863 
1864  v += sigma*det;
1865  vcell += sigma*det;
1866  }
1867  }
1868  }
1869  if (vmin > vcell)
1870  vmin = vcell;
1871  }
1872  }
1873  }
1874 
1875  // Fill in return value references
1876  vol = v/static_cast<double>(_n_cells);
1877  Vmin = vmin;
1878 
1879  return gqmin;
1880 }
1881 
1882 
1883 
1884 // Executes one step of minimization algorithm:
1885 // finds minimization direction (P=H^{-1} \grad J) and solves approximately
1886 // local minimization problem for optimal step in this minimization direction (tau=min J(R+tau P))
1887 double VariationalMeshSmoother::minJ(Array2D<double> & R,
1888  const std::vector<int> & mask,
1889  const Array2D<int> & cells,
1890  const std::vector<int> & mcells,
1891  double epsilon,
1892  double w,
1893  int me,
1894  const Array3D<double> & H,
1895  double vol,
1896  const std::vector<int> & edges,
1897  const std::vector<int> & hnodes,
1898  int msglev,
1899  double & Vmin,
1900  double & emax,
1901  double & qmin,
1902  int adp,
1903  const std::vector<double> & afun)
1904 {
1905  // columns - max number of nonzero entries in every row of global matrix
1906  int columns = _dim*_dim*10;
1907 
1908  // local Hessian matrix
1909  Array3D<double> W(_dim, 3*_dim + _dim%2, 3*_dim + _dim%2);
1910 
1911  // F - local gradient
1912  Array2D<double> F(_dim, 3*_dim + _dim%2);
1913 
1914  Array2D<double> Rpr(_n_nodes, _dim);
1915 
1916  // P - minimization direction
1917  Array2D<double> P(_n_nodes, _dim);
1918 
1919  // A, JA - internal form of global matrix
1920  Array2D<int> JA(_dim*_n_nodes, columns);
1921  Array2D<double> A(_dim*_n_nodes, columns);
1922 
1923  // G - adaptation metric
1924  Array2D<double> G(_n_cells, _dim);
1925 
1926  // rhs for solver
1927  std::vector<double> b(_dim*_n_nodes);
1928 
1929  // u - solution vector
1930  std::vector<double> u(_dim*_n_nodes);
1931 
1932  // matrix
1933  std::vector<double> a(_dim*_n_nodes*columns);
1934  std::vector<int> ia(_dim*_n_nodes + 1);
1935  std::vector<int> ja(_dim*_n_nodes*columns);
1936 
1937  // nonzero - norm of gradient
1938  double nonzero = 0.;
1939 
1940  // Jpr - value of functional
1941  double Jpr = 0.;
1942 
1943  // find minimization direction P
1944  for (dof_id_type i=0; i<_n_cells; i++)
1945  {
1946  int nvert = 0;
1947  while (cells[i][nvert] >= 0)
1948  nvert++;
1949 
1950  // determination of local matrices on each cell
1951  for (unsigned j=0; j<_dim; j++)
1952  {
1953  G[i][j] = 0; // adaptation metric G is held constant throughout minJ run
1954  if (adp < 0)
1955  {
1956  for (int k=0; k<std::abs(adp); k++)
1957  G[i][j] += afun[i*(-adp)+k]; // cell-based adaptivity is computed here
1958  }
1959  }
1960  for (unsigned index=0; index<_dim; index++)
1961  {
1962  // initialize local matrices
1963  for (unsigned k=0; k<3*_dim + _dim%2; k++)
1964  {
1965  F[index][k] = 0;
1966 
1967  for (unsigned j=0; j<3*_dim + _dim%2; j++)
1968  W[index][k][j] = 0;
1969  }
1970  }
1971  if (mcells[i] >= 0)
1972  {
1973  // if cell is not excluded
1974  double lVmin, lqmin;
1975  Jpr += localP(W, F, R, cells[i], mask, epsilon, w, nvert, H[i],
1976  me, vol, 0, lVmin, lqmin, adp, afun, G[i]);
1977  }
1978  else
1979  {
1980  for (unsigned index=0; index<_dim; index++)
1981  for (int j=0; j<nvert; j++)
1982  W[index][j][j] = 1;
1983  }
1984 
1985  // assembly of an upper triangular part of a global matrix A
1986  for (unsigned index=0; index<_dim; index++)
1987  {
1988  for (int l=0; l<nvert; l++)
1989  {
1990  for (int m=0; m<nvert; m++)
1991  {
1992  if ((W[index][l][m] != 0) &&
1993  (cells[i][m] >= cells[i][l]))
1994  {
1995  int sch = 0;
1996  int ind = 1;
1997  while (ind != 0)
1998  {
1999  if (A[cells[i][l] + index*_n_nodes][sch] != 0)
2000  {
2001  if (JA[cells[i][l] + index*_n_nodes][sch] == static_cast<int>(cells[i][m] + index*_n_nodes))
2002  {
2003  A[cells[i][l] + index*_n_nodes][sch] = A[cells[i][l] + index*_n_nodes][sch] + W[index][l][m];
2004  ind=0;
2005  }
2006  else
2007  sch++;
2008  }
2009  else
2010  {
2011  A[cells[i][l] + index*_n_nodes][sch] = W[index][l][m];
2012  JA[cells[i][l] + index*_n_nodes][sch] = cells[i][m] + index*_n_nodes;
2013  ind = 0;
2014  }
2015 
2016  if (sch > columns-1)
2017  _logfile << "error: # of nonzero entries in the "
2018  << cells[i][l]
2019  << " row of Hessian ="
2020  << sch
2021  << ">= columns="
2022  << columns
2023  << std::endl;
2024  }
2025  }
2026  }
2027  b[cells[i][l] + index*_n_nodes] = b[cells[i][l] + index*_n_nodes] - F[index][l];
2028  }
2029  }
2030  // end of matrix A
2031  }
2032 
2033  // HN correction
2034 
2035  // tolerance for HN being the mid-edge node for its parents
2036  double Tau_hn = pow(vol, 1./static_cast<double>(_dim))*1e-10;
2037  for (dof_id_type i=0; i<_n_hanging_edges; i++)
2038  {
2039  int ind_i = hnodes[i];
2040  int ind_j = edges[2*i];
2041  int ind_k = edges[2*i+1];
2042 
2043  for (unsigned j=0; j<_dim; j++)
2044  {
2045  int g_i = int(R[ind_i][j] - 0.5*(R[ind_j][j]+R[ind_k][j]));
2046  Jpr += g_i*g_i/(2*Tau_hn);
2047  b[ind_i + j*_n_nodes] -= g_i/Tau_hn;
2048  b[ind_j + j*_n_nodes] += 0.5*g_i/Tau_hn;
2049  b[ind_k + j*_n_nodes] += 0.5*g_i/Tau_hn;
2050  }
2051 
2052  for (int ind=0; ind<columns; ind++)
2053  {
2054  if (JA[ind_i][ind] == ind_i)
2055  A[ind_i][ind] += 1./Tau_hn;
2056 
2057  if (JA[ind_i+_n_nodes][ind] == static_cast<int>(ind_i+_n_nodes))
2058  A[ind_i+_n_nodes][ind] += 1./Tau_hn;
2059 
2060  if (_dim == 3)
2061  if (JA[ind_i+2*_n_nodes][ind] == static_cast<int>(ind_i+2*_n_nodes))
2062  A[ind_i+2*_n_nodes][ind] += 1./Tau_hn;
2063 
2064  if ((JA[ind_i][ind] == ind_j) ||
2065  (JA[ind_i][ind] == ind_k))
2066  A[ind_i][ind] -= 0.5/Tau_hn;
2067 
2068  if ((JA[ind_i+_n_nodes][ind] == static_cast<int>(ind_j+_n_nodes)) ||
2069  (JA[ind_i+_n_nodes][ind] == static_cast<int>(ind_k+_n_nodes)))
2070  A[ind_i+_n_nodes][ind] -= 0.5/Tau_hn;
2071 
2072  if (_dim == 3)
2073  if ((JA[ind_i+2*_n_nodes][ind] == static_cast<int>(ind_j+2*_n_nodes)) ||
2074  (JA[ind_i+2*_n_nodes][ind] == static_cast<int>(ind_k+2*_n_nodes)))
2075  A[ind_i+2*_n_nodes][ind] -= 0.5/Tau_hn;
2076 
2077  if (JA[ind_j][ind] == ind_i)
2078  A[ind_j][ind] -= 0.5/Tau_hn;
2079 
2080  if (JA[ind_k][ind] == ind_i)
2081  A[ind_k][ind] -= 0.5/Tau_hn;
2082 
2083  if (JA[ind_j+_n_nodes][ind] == static_cast<int>(ind_i+_n_nodes))
2084  A[ind_j+_n_nodes][ind] -= 0.5/Tau_hn;
2085 
2086  if (JA[ind_k+_n_nodes][ind] == static_cast<int>(ind_i+_n_nodes))
2087  A[ind_k+_n_nodes][ind] -= 0.5/Tau_hn;
2088 
2089  if (_dim == 3)
2090  if (JA[ind_j+2*_n_nodes][ind] == static_cast<int>(ind_i+2*_n_nodes))
2091  A[ind_j+2*_n_nodes][ind] -= 0.5/Tau_hn;
2092 
2093  if (_dim == 3)
2094  if (JA[ind_k+2*_n_nodes][ind] == static_cast<int>(ind_i+2*_n_nodes))
2095  A[ind_k+2*_n_nodes][ind] -= 0.5/Tau_hn;
2096 
2097  if ((JA[ind_j][ind] == ind_j) ||
2098  (JA[ind_j][ind] == ind_k))
2099  A[ind_j][ind] += 0.25/Tau_hn;
2100 
2101  if ((JA[ind_k][ind] == ind_j) ||
2102  (JA[ind_k][ind] == ind_k))
2103  A[ind_k][ind] += 0.25/Tau_hn;
2104 
2105  if ((JA[ind_j+_n_nodes][ind] == static_cast<int>(ind_j+_n_nodes)) ||
2106  (JA[ind_j+_n_nodes][ind] == static_cast<int>(ind_k+_n_nodes)))
2107  A[ind_j+_n_nodes][ind] += 0.25/Tau_hn;
2108 
2109  if ((JA[ind_k+_n_nodes][ind] == static_cast<int>(ind_j+_n_nodes)) ||
2110  (JA[ind_k+_n_nodes][ind] == static_cast<int>(ind_k+_n_nodes)))
2111  A[ind_k+_n_nodes][ind] += 0.25/Tau_hn;
2112 
2113  if (_dim == 3)
2114  if ((JA[ind_j+2*_n_nodes][ind] == static_cast<int>(ind_j+2*_n_nodes)) ||
2115  (JA[ind_j+2*_n_nodes][ind] == static_cast<int>(ind_k+2*_n_nodes)))
2116  A[ind_j+2*_n_nodes][ind] += 0.25/Tau_hn;
2117 
2118  if (_dim==3)
2119  if ((JA[ind_k+2*_n_nodes][ind] == static_cast<int>(ind_j+2*_n_nodes)) ||
2120  (JA[ind_k+2*_n_nodes][ind] == static_cast<int>(ind_k+2*_n_nodes)))
2121  A[ind_k+2*_n_nodes][ind] += 0.25/Tau_hn;
2122  }
2123  }
2124 
2125  // ||\grad J||_2
2126  for (dof_id_type i=0; i<_dim*_n_nodes; i++)
2127  nonzero += b[i]*b[i];
2128 
2129  // sort matrix A
2130  for (dof_id_type i=0; i<_dim*_n_nodes; i++)
2131  {
2132  for (int j=columns-1; j>1; j--)
2133  {
2134  for (int k=0; k<j; k++)
2135  {
2136  if (JA[i][k] > JA[i][k+1])
2137  {
2138  int sch = JA[i][k];
2139  JA[i][k] = JA[i][k+1];
2140  JA[i][k+1] = sch;
2141  double tmp = A[i][k];
2142  A[i][k] = A[i][k+1];
2143  A[i][k+1] = tmp;
2144  }
2145  }
2146  }
2147  }
2148 
2149  double eps = std::sqrt(vol)*1e-9;
2150 
2151  // solver for P (unconstrained)
2152  ia[0] = 0;
2153  {
2154  int j = 0;
2155  for (dof_id_type i=0; i<_dim*_n_nodes; i++)
2156  {
2157  u[i] = 0;
2158  int nz = 0;
2159  for (int k=0; k<columns; k++)
2160  {
2161  if (A[i][k] != 0)
2162  {
2163  nz++;
2164  ja[j] = JA[i][k]+1;
2165  a[j] = A[i][k];
2166  j++;
2167  }
2168  }
2169  ia[i+1] = ia[i] + nz;
2170  }
2171  }
2172 
2173  dof_id_type m = _dim*_n_nodes;
2174  int sch = (msglev >= 3) ? 1 : 0;
2175 
2176  solver(m, ia, ja, a, u, b, eps, 100, sch);
2177  // sol_pcg_pj(m, ia, ja, a, u, b, eps, 100, sch);
2178 
2179  for (dof_id_type i=0; i<_n_nodes; i++)
2180  {
2181  //ensure fixed nodes are not moved
2182  for (unsigned index=0; index<_dim; index++)
2183  if (mask[i] == 1)
2184  P[i][index] = 0;
2185  else
2186  P[i][index] = u[i+index*_n_nodes];
2187  }
2188 
2189  // P is determined
2190  if (msglev >= 4)
2191  {
2192  for (dof_id_type i=0; i<_n_nodes; i++)
2193  {
2194  for (unsigned j=0; j<_dim; j++)
2195  if (P[i][j] != 0)
2196  _logfile << "P[" << i << "][" << j << "]=" << P[i][j];
2197  }
2198  }
2199 
2200  // local minimization problem, determination of tau
2201  if (msglev >= 3)
2202  _logfile << "dJ=" << std::sqrt(nonzero) << " J0=" << Jpr << std::endl;
2203 
2204  double
2205  J = 1.e32,
2206  tau = 0.,
2207  gVmin = 0.,
2208  gemax = 0.,
2209  gqmin = 0.;
2210 
2211  int j = 1;
2212 
2213  while ((Jpr <= J) && (j > -30))
2214  {
2215  j = j-1;
2216 
2217  // search through finite # of values for tau
2218  tau = pow(2., j);
2219  for (dof_id_type i=0; i<_n_nodes; i++)
2220  for (unsigned k=0; k<_dim; k++)
2221  Rpr[i][k] = R[i][k] + tau*P[i][k];
2222 
2223  J = 0;
2224  gVmin = 1e32;
2225  gemax = -1e32;
2226  gqmin = 1e32;
2227  for (dof_id_type i=0; i<_n_cells; i++)
2228  {
2229  if (mcells[i] >= 0)
2230  {
2231  int nvert = 0;
2232  while (cells[i][nvert] >= 0)
2233  nvert++;
2234 
2235  double lVmin, lqmin;
2236  double lemax = localP(W, F, Rpr, cells[i], mask, epsilon, w, nvert, H[i], me, vol, 1, lVmin, lqmin, adp, afun, G[i]);
2237 
2238  J += lemax;
2239  if (gVmin > lVmin)
2240  gVmin = lVmin;
2241 
2242  if (gemax < lemax)
2243  gemax = lemax;
2244 
2245  if (gqmin > lqmin)
2246  gqmin = lqmin;
2247 
2248  // HN correction
2249  for (dof_id_type ii=0; ii<_n_hanging_edges; ii++)
2250  {
2251  int ind_i = hnodes[ii];
2252  int ind_j = edges[2*ii];
2253  int ind_k = edges[2*ii+1];
2254  for (unsigned jj=0; jj<_dim; jj++)
2255  {
2256  int g_i = int(Rpr[ind_i][jj] - 0.5*(Rpr[ind_j][jj]+Rpr[ind_k][jj]));
2257  J += g_i*g_i/(2*Tau_hn);
2258  }
2259  }
2260  }
2261  }
2262  if (msglev >= 3)
2263  _logfile << "tau=" << tau << " J=" << J << std::endl;
2264  }
2265 
2266 
2267  double
2268  T = 0.,
2269  gtmin0 = 0.,
2270  gtmax0 = 0.,
2271  gqmin0 = 0.;
2272 
2273  if (j != -30)
2274  {
2275  Jpr = J;
2276  for (unsigned i=0; i<_n_nodes; i++)
2277  for (unsigned k=0; k<_dim; k++)
2278  Rpr[i][k] = R[i][k] + tau*0.5*P[i][k];
2279 
2280  J = 0;
2281  gtmin0 = 1e32;
2282  gtmax0 = -1e32;
2283  gqmin0 = 1e32;
2284  for (dof_id_type i=0; i<_n_cells; i++)
2285  {
2286  if (mcells[i] >= 0)
2287  {
2288  int nvert = 0;
2289  while (cells[i][nvert] >= 0)
2290  nvert++;
2291 
2292  double lVmin, lqmin;
2293  double lemax = localP(W, F, Rpr, cells[i], mask, epsilon, w, nvert, H[i], me, vol, 1, lVmin, lqmin, adp, afun, G[i]);
2294  J += lemax;
2295 
2296  if (gtmin0 > lVmin)
2297  gtmin0 = lVmin;
2298 
2299  if (gtmax0 < lemax)
2300  gtmax0 = lemax;
2301 
2302  if (gqmin0 > lqmin)
2303  gqmin0 = lqmin;
2304 
2305  // HN correction
2306  for (dof_id_type ii=0; ii<_n_hanging_edges; ii++)
2307  {
2308  int ind_i = hnodes[ii];
2309  int ind_j = edges[2*ii];
2310  int ind_k = edges[2*ii+1];
2311 
2312  for (unsigned jj=0; jj<_dim; jj++)
2313  {
2314  int g_i = int(Rpr[ind_i][jj] - 0.5*(Rpr[ind_j][jj] + Rpr[ind_k][jj]));
2315  J += g_i*g_i/(2*Tau_hn);
2316  }
2317  }
2318  }
2319  }
2320  }
2321 
2322  if (Jpr > J)
2323  {
2324  T = 0.5*tau;
2325  // Set up return values passed by reference
2326  Vmin = gtmin0;
2327  emax = gtmax0;
2328  qmin = gqmin0;
2329  }
2330  else
2331  {
2332  T = tau;
2333  J = Jpr;
2334  // Set up return values passed by reference
2335  Vmin = gVmin;
2336  emax = gemax;
2337  qmin = gqmin;
2338  }
2339 
2340  nonzero = 0;
2341  for (dof_id_type j2=0; j2<_n_nodes; j2++)
2342  for (unsigned k=0; k<_dim; k++)
2343  {
2344  R[j2][k] = R[j2][k] + T*P[j2][k];
2345  nonzero += T*P[j2][k]*T*P[j2][k];
2346  }
2347 
2348  if (msglev >= 2)
2349  _logfile << "tau=" << T << ", J=" << J << std::endl;
2350 
2351  return std::sqrt(nonzero);
2352 }
2353 
2354 
2355 
2356 // minJ() with sliding Boundary Nodes constraints and no account for HN,
2357 // using Lagrange multiplier formulation: minimize L=J+\sum lam*g;
2358 // only works in 2D
2359 double VariationalMeshSmoother::minJ_BC(Array2D<double> & R,
2360  const std::vector<int> & mask,
2361  const Array2D<int> & cells,
2362  const std::vector<int> & mcells,
2363  double epsilon,
2364  double w,
2365  int me,
2366  const Array3D<double> & H,
2367  double vol,
2368  int msglev,
2369  double & Vmin,
2370  double & emax,
2371  double & qmin,
2372  int adp,
2373  const std::vector<double> & afun,
2374  int NCN)
2375 {
2376  // new form of matrices, 5 iterations for minL
2377  double tau = 0., J = 0., T, Jpr, L, lVmin, lqmin, gVmin = 0., gqmin = 0., gVmin0 = 0.,
2378  gqmin0 = 0., lemax, gemax = 0., gemax0 = 0.;
2379 
2380  // array of sliding BN
2381  std::vector<int> Bind(NCN);
2382  std::vector<double> lam(NCN);
2383  std::vector<double> hm(2*_n_nodes);
2384  std::vector<double> Plam(NCN);
2385 
2386  // holds constraints = local approximation to the boundary
2387  std::vector<double> constr(4*NCN);
2388 
2389  Array2D<double> F(2, 6);
2390  Array3D<double> W(2, 6, 6);
2391  Array2D<double> Rpr(_n_nodes, 2);
2392  Array2D<double> P(_n_nodes, 2);
2393 
2394  std::vector<double> b(2*_n_nodes);
2395 
2396  Array2D<double> G(_n_cells, 6);
2397 
2398  // assembler of constraints
2399  const double eps = std::sqrt(vol)*1e-9;
2400 
2401  for (int i=0; i<4*NCN; i++)
2402  constr[i] = 1./eps;
2403 
2404  {
2405  int I = 0;
2406  for (dof_id_type i=0; i<_n_nodes; i++)
2407  if (mask[i] == 2)
2408  {
2409  Bind[I] = i;
2410  I++;
2411  }
2412  }
2413 
2414  for (int I=0; I<NCN; I++)
2415  {
2416  // The boundary connectivity loop sets up the j and k indices
2417  // but I am not sure about the logic of this: j and k are overwritten
2418  // at every iteration of the boundary connectivity loop and only used
2419  // *after* the boundary connectivity loop -- this seems like a bug but
2420  // I maintained the original behavior of the algorithm [JWP].
2421  int
2422  i = Bind[I],
2423  j = 0,
2424  k = 0,
2425  ind = 0;
2426 
2427  // boundary connectivity
2428  for (dof_id_type l=0; l<_n_cells; l++)
2429  {
2430  int nvert = 0;
2431 
2432  while (cells[l][nvert] >= 0)
2433  nvert++;
2434 
2435  switch (nvert)
2436  {
2437  case 3: // tri
2438  if (i == cells[l][0])
2439  {
2440  if (mask[cells[l][1]] > 0)
2441  {
2442  if (ind == 0)
2443  j = cells[l][1];
2444  else
2445  k = cells[l][1];
2446  ind++;
2447  }
2448 
2449  if (mask[cells[l][2]] > 0)
2450  {
2451  if (ind == 0)
2452  j = cells[l][2];
2453  else
2454  k = cells[l][2];
2455  ind++;
2456  }
2457  }
2458 
2459  if (i == cells[l][1])
2460  {
2461  if (mask[cells[l][0]] > 0)
2462  {
2463  if (ind == 0)
2464  j = cells[l][0];
2465  else
2466  k = cells[l][0];
2467  ind++;
2468  }
2469 
2470  if (mask[cells[l][2]] > 0)
2471  {
2472  if (ind == 0)
2473  j = cells[l][2];
2474  else
2475  k = cells[l][2];
2476  ind++;
2477  }
2478  }
2479 
2480  if (i == cells[l][2])
2481  {
2482  if (mask[cells[l][1]] > 0)
2483  {
2484  if (ind == 0)
2485  j = cells[l][1];
2486  else
2487  k = cells[l][1];
2488  ind++;
2489  }
2490 
2491  if (mask[cells[l][0]] > 0)
2492  {
2493  if (ind == 0)
2494  j = cells[l][0];
2495  else
2496  k = cells[l][0];
2497  ind++;
2498  }
2499  }
2500  break;
2501 
2502  case 4: // quad
2503  if ((i == cells[l][0]) ||
2504  (i == cells[l][3]))
2505  {
2506  if (mask[cells[l][1]] > 0)
2507  {
2508  if (ind == 0)
2509  j = cells[l][1];
2510  else
2511  k = cells[l][1];
2512  ind++;
2513  }
2514 
2515  if (mask[cells[l][2]] > 0)
2516  {
2517  if (ind == 0)
2518  j = cells[l][2];
2519  else
2520  k = cells[l][2];
2521  ind++;
2522  }
2523  }
2524 
2525  if ((i == cells[l][1]) ||
2526  (i == cells[l][2]))
2527  {
2528  if (mask[cells[l][0]] > 0)
2529  {
2530  if (ind == 0)
2531  j = cells[l][0];
2532  else
2533  k = cells[l][0];
2534  ind++;
2535  }
2536 
2537  if (mask[cells[l][3]] > 0)
2538  {
2539  if (ind == 0)
2540  j = cells[l][3];
2541  else
2542  k = cells[l][3];
2543  ind++;
2544  }
2545  }
2546  break;
2547 
2548  case 6: //quad tri
2549  if (i == cells[l][0])
2550  {
2551  if (mask[cells[l][1]] > 0)
2552  {
2553  if (ind == 0)
2554  j = cells[l][5];
2555  else
2556  k = cells[l][5];
2557  ind++;
2558  }
2559 
2560  if (mask[cells[l][2]] > 0)
2561  {
2562  if (ind == 0)
2563  j = cells[l][4];
2564  else
2565  k = cells[l][4];
2566  ind++;
2567  }
2568  }
2569 
2570  if (i == cells[l][1])
2571  {
2572  if (mask[cells[l][0]] > 0)
2573  {
2574  if (ind == 0)
2575  j = cells[l][5];
2576  else
2577  k = cells[l][5];
2578  ind++;
2579  }
2580 
2581  if (mask[cells[l][2]] > 0)
2582  {
2583  if (ind == 0)
2584  j = cells[l][3];
2585  else
2586  k = cells[l][3];
2587  ind++;
2588  }
2589  }
2590 
2591  if (i == cells[l][2])
2592  {
2593  if (mask[cells[l][1]] > 0)
2594  {
2595  if (ind == 0)
2596  j = cells[l][3];
2597  else
2598  k = cells[l][3];
2599  ind++;
2600  }
2601 
2602  if (mask[cells[l][0]] > 0)
2603  {
2604  if (ind == 0)
2605  j = cells[l][4];
2606  else
2607  k = cells[l][4];
2608  ind++;
2609  }
2610  }
2611 
2612  if (i == cells[l][3])
2613  {
2614  j = cells[l][1];
2615  k = cells[l][2];
2616  }
2617 
2618  if (i == cells[l][4])
2619  {
2620  j = cells[l][2];
2621  k = cells[l][0];
2622  }
2623 
2624  if (i == cells[l][5])
2625  {
2626  j = cells[l][0];
2627  k = cells[l][1];
2628  }
2629  break;
2630 
2631  default:
2632  libmesh_error_msg("Unrecognized nvert = " << nvert);
2633  }
2634  } // end boundary connectivity
2635 
2636  // lines
2637  double del1 = R[j][0] - R[i][0];
2638  double del2 = R[i][0] - R[k][0];
2639 
2640  if ((std::abs(del1) < eps) &&
2641  (std::abs(del2) < eps))
2642  {
2643  constr[I*4] = 1;
2644  constr[I*4+1] = 0;
2645  constr[I*4+2] = R[i][0];
2646  constr[I*4+3] = R[i][1];
2647  }
2648  else
2649  {
2650  del1 = R[j][1] - R[i][1];
2651  del2 = R[i][1] - R[k][1];
2652  if ((std::abs(del1) < eps) &&
2653  (std::abs(del2) < eps))
2654  {
2655  constr[I*4] = 0;
2656  constr[I*4+1] = 1;
2657  constr[I*4+2] = R[i][0];
2658  constr[I*4+3] = R[i][1];
2659  }
2660  else
2661  {
2662  J =
2663  (R[i][0]-R[j][0])*(R[k][1]-R[j][1]) -
2664  (R[k][0]-R[j][0])*(R[i][1]-R[j][1]);
2665 
2666  if (std::abs(J) < eps)
2667  {
2668  constr[I*4] = 1./(R[k][0]-R[j][0]);
2669  constr[I*4+1] = -1./(R[k][1]-R[j][1]);
2670  constr[I*4+2] = R[i][0];
2671  constr[I*4+3] = R[i][1];
2672  }
2673  else
2674  {
2675  // circle
2676  double x0 = ((R[k][1]-R[j][1])*(R[i][0]*R[i][0]-R[j][0]*R[j][0]+R[i][1]*R[i][1]-R[j][1]*R[j][1]) +
2677  (R[j][1]-R[i][1])*(R[k][0]*R[k][0]-R[j][0]*R[j][0]+R[k][1]*R[k][1]-R[j][1]*R[j][1]))/(2*J);
2678  double y0 = ((R[j][0]-R[k][0])*(R[i][0]*R[i][0]-R[j][0]*R[j][0]+R[i][1]*R[i][1]-R[j][1]*R[j][1]) +
2679  (R[i][0]-R[j][0])*(R[k][0]*R[k][0]-R[j][0]*R[j][0]+R[k][1]*R[k][1]-R[j][1]*R[j][1]))/(2*J);
2680  constr[I*4] = x0;
2681  constr[I*4+1] = y0;
2682  constr[I*4+2] = (R[i][0]-x0)*(R[i][0]-x0) + (R[i][1]-y0)*(R[i][1]-y0);
2683  }
2684  }
2685  }
2686  }
2687 
2688  // for (int i=0; i<NCN; i++){
2689  // for (int j=0; j<4; j++) fprintf(stdout," %e ",constr[4*i+j]);
2690  // fprintf(stdout, "constr %d node %d \n",i,Bind[i]);
2691  // }
2692 
2693  // initial values
2694  for (int i=0; i<NCN; i++)
2695  lam[i] = 0;
2696 
2697  // Eventual return value
2698  double nonzero = 0.;
2699 
2700  // Temporary result variable
2701  double qq = 0.;
2702 
2703  for (int nz=0; nz<5; nz++)
2704  {
2705  // find H and -grad J
2706  nonzero = 0.;
2707  Jpr = 0;
2708  for (dof_id_type i=0; i<2*_n_nodes; i++)
2709  {
2710  b[i] = 0;
2711  hm[i] = 0;
2712  }
2713 
2714  for (dof_id_type i=0; i<_n_cells; i++)
2715  {
2716  int nvert = 0;
2717 
2718  while (cells[i][nvert] >= 0)
2719  nvert++;
2720 
2721  for (int j=0; j<nvert; j++)
2722  {
2723  G[i][j] = 0;
2724  if (adp < 0)
2725  for (int k=0; k<std::abs(adp); k++)
2726  G[i][j] += afun[i*(-adp) + k];
2727  }
2728 
2729  for (int index=0; index<2; index++)
2730  for (int k=0; k<nvert; k++)
2731  {
2732  F[index][k] = 0;
2733  for (int j=0; j<nvert; j++)
2734  W[index][k][j] = 0;
2735  }
2736 
2737  if (mcells[i] >= 0)
2738  Jpr += localP(W, F, R, cells[i], mask, epsilon, w, nvert, H[i],
2739  me, vol, 0, lVmin, lqmin, adp, afun, G[i]);
2740 
2741  else
2742  {
2743  for (unsigned index=0; index<2; index++)
2744  for (int j=0; j<nvert; j++)
2745  W[index][j][j] = 1;
2746  }
2747 
2748  for (unsigned index=0; index<2; index++)
2749  for (int l=0; l<nvert; l++)
2750  {
2751  // diagonal Hessian
2752  hm[cells[i][l] + index*_n_nodes] += W[index][l][l];
2753  b[cells[i][l] + index*_n_nodes] -= F[index][l];
2754  }
2755  }
2756 
2757  // ||grad J||_2
2758  for (dof_id_type i=0; i<2*_n_nodes; i++)
2759  nonzero += b[i]*b[i];
2760 
2761  // solve for Plam
2762  for (int I=0; I<NCN; I++)
2763  {
2764  int i = Bind[I];
2765  double
2766  Bx = 0.,
2767  By = 0.,
2768  g = 0.;
2769 
2770  if (constr[4*I+3] < 0.5/eps)
2771  {
2772  Bx = constr[4*I];
2773  By = constr[4*I+1];
2774  g = (R[i][0]-constr[4*I+2])*constr[4*I] + (R[i][1]-constr[4*I+3])*constr[4*I+1];
2775  }
2776  else
2777  {
2778  Bx = 2*(R[i][0] - constr[4*I]);
2779  By = 2*(R[i][1] - constr[4*I+1]);
2780  hm[i] += 2*lam[I];
2781  hm[i+_n_nodes] += 2*lam[I];
2782  g =
2783  (R[i][0] - constr[4*I])*(R[i][0]-constr[4*I]) +
2784  (R[i][1] - constr[4*I+1])*(R[i][1]-constr[4*I+1]) - constr[4*I+2];
2785  }
2786 
2787  Jpr += lam[I]*g;
2788  qq = Bx*b[i]/hm[i] + By*b[i+_n_nodes]/hm[i+_n_nodes] - g;
2789  double a = Bx*Bx/hm[i] + By*By/hm[i+_n_nodes];
2790 
2791  if (a != 0)
2792  Plam[I] = qq/a;
2793  else
2794  _logfile << "error: B^TH-1B is degenerate" << std::endl;
2795 
2796  b[i] -= Plam[I]*Bx;
2797  b[i+_n_nodes] -= Plam[I]*By;
2798  Plam[I] -= lam[I];
2799  }
2800 
2801  // solve for P
2802  for (dof_id_type i=0; i<_n_nodes; i++)
2803  {
2804  P[i][0] = b[i]/hm[i];
2805  P[i][1] = b[i+_n_nodes]/hm[i+_n_nodes];
2806  }
2807 
2808  // correct solution
2809  for (dof_id_type i=0; i<_n_nodes; i++)
2810  for (unsigned j=0; j<2; j++)
2811  if ((std::abs(P[i][j]) < eps) || (mask[i] == 1))
2812  P[i][j] = 0;
2813 
2814  // P is determined
2815  if (msglev >= 3)
2816  {
2817  for (dof_id_type i=0; i<_n_nodes; i++)
2818  for (unsigned j=0; j<2; j++)
2819  if (P[i][j] != 0)
2820  _logfile << "P[" << i << "][" << j << "]=" << P[i][j] << " ";
2821  }
2822 
2823  // local minimization problem, determination of tau
2824  if (msglev >= 3)
2825  _logfile << "dJ=" << std::sqrt(nonzero) << " L0=" << Jpr << std::endl;
2826 
2827  L = 1.e32;
2828  int j = 1;
2829 
2830  while ((Jpr <= L) && (j > -30))
2831  {
2832  j = j-1;
2833  tau = pow(2.,j);
2834  for (dof_id_type i=0; i<_n_nodes; i++)
2835  for (unsigned k=0; k<2; k++)
2836  Rpr[i][k] = R[i][k] + tau*P[i][k];
2837 
2838  J = 0;
2839  gVmin = 1.e32;
2840  gemax = -1.e32;
2841  gqmin = 1.e32;
2842  for (dof_id_type i=0; i<_n_cells; i++)
2843  if (mcells[i] >= 0)
2844  {
2845  int nvert = 0;
2846  while (cells[i][nvert] >= 0)
2847  nvert++;
2848 
2849  lemax = localP(W, F, Rpr, cells[i], mask, epsilon, w, nvert, H[i], me, vol, 1, lVmin,
2850  lqmin, adp, afun, G[i]);
2851  J += lemax;
2852 
2853  if (gVmin > lVmin)
2854  gVmin = lVmin;
2855 
2856  if (gemax < lemax)
2857  gemax = lemax;
2858 
2859  if (gqmin > lqmin)
2860  gqmin = lqmin;
2861  }
2862 
2863  L = J;
2864 
2865  // constraints contribution
2866  for (int I=0; I<NCN; I++)
2867  {
2868  int i = Bind[I];
2869  double g = 0.;
2870 
2871  if (constr[4*I+3] < 0.5/eps)
2872  g = (Rpr[i][0] - constr[4*I+2])*constr[4*I] + (Rpr[i][1]-constr[4*I+3])*constr[4*I+1];
2873 
2874  else
2875  g = (Rpr[i][0]-constr[4*I])*(Rpr[i][0]-constr[4*I]) +
2876  (Rpr[i][1]-constr[4*I+1])*(Rpr[i][1]-constr[4*I+1]) - constr[4*I+2];
2877 
2878  L += (lam[I] + tau*Plam[I])*g;
2879  }
2880 
2881  // end of constraints
2882  if (msglev >= 3)
2883  _logfile << " tau=" << tau << " J=" << J << std::endl;
2884  } // end while
2885 
2886  if (j == -30)
2887  T = 0;
2888  else
2889  {
2890  Jpr = L;
2891  qq = J;
2892  for (dof_id_type i=0; i<_n_nodes; i++)
2893  for (unsigned k=0; k<2; k++)
2894  Rpr[i][k] = R[i][k] + tau*0.5*P[i][k];
2895 
2896  J = 0;
2897  gVmin0 = 1.e32;
2898  gemax0 = -1.e32;
2899  gqmin0 = 1.e32;
2900 
2901  for (dof_id_type i=0; i<_n_cells; i++)
2902  if (mcells[i] >= 0)
2903  {
2904  int nvert = 0;
2905  while (cells[i][nvert] >= 0)
2906  nvert++;
2907 
2908  lemax = localP(W, F, Rpr, cells[i], mask, epsilon, w, nvert, H[i], me, vol, 1, lVmin,
2909  lqmin, adp, afun, G[i]);
2910  J += lemax;
2911 
2912  if (gVmin0 > lVmin)
2913  gVmin0 = lVmin;
2914 
2915  if (gemax0 < lemax)
2916  gemax0 = lemax;
2917 
2918  if (gqmin0 > lqmin)
2919  gqmin0 = lqmin;
2920  }
2921 
2922  L = J;
2923 
2924  // constraints contribution
2925  for (int I=0; I<NCN; I++)
2926  {
2927  int i = Bind[I];
2928  double g = 0.;
2929 
2930  if (constr[4*I+3] < 0.5/eps)
2931  g = (Rpr[i][0]-constr[4*I+2])*constr[4*I] + (Rpr[i][1]-constr[4*I+3])*constr[4*I+1];
2932 
2933  else
2934  g = (Rpr[i][0]-constr[4*I])*(Rpr[i][0]-constr[4*I]) +
2935  (Rpr[i][1]-constr[4*I+1])*(Rpr[i][1]-constr[4*I+1]) - constr[4*I+2];
2936 
2937  L += (lam[I] + tau*0.5*Plam[I])*g;
2938  }
2939  // end of constraints
2940  }
2941 
2942  if (Jpr > L)
2943  {
2944  T = 0.5*tau;
2945  // Set return values by reference
2946  Vmin = gVmin0;
2947  emax = gemax0;
2948  qmin = gqmin0;
2949  }
2950  else
2951  {
2952  T = tau;
2953  L = Jpr;
2954  J = qq;
2955  // Set return values by reference
2956  Vmin = gVmin;
2957  emax = gemax;
2958  qmin = gqmin;
2959  }
2960 
2961  for (dof_id_type i=0; i<_n_nodes; i++)
2962  for (unsigned k=0; k<2; k++)
2963  R[i][k] += T*P[i][k];
2964 
2965  for (int i=0; i<NCN; i++)
2966  lam[i] += T*Plam[i];
2967 
2968  } // end Lagrangian iter
2969 
2970  if (msglev >= 2)
2971  _logfile << "tau=" << T << ", J=" << J << ", L=" << L << std::endl;
2972 
2973  return std::sqrt(nonzero);
2974 }
2975 
2976 
2977 
2978 // composes local matrix W and right side F from all quadrature nodes of one cell
2979 double VariationalMeshSmoother::localP(Array3D<double> & W,
2980  Array2D<double> & F,
2981  Array2D<double> & R,
2982  const std::vector<int> & cell_in,
2983  const std::vector<int> & mask,
2984  double epsilon,
2985  double w,
2986  int nvert,
2987  const Array2D<double> & H,
2988  int me,
2989  double vol,
2990  int f,
2991  double & Vmin,
2992  double & qmin,
2993  int adp,
2994  const std::vector<double> & afun,
2995  std::vector<double> & Gloc)
2996 {
2997  // K - determines approximation rule for local integral over the cell
2998  std::vector<double> K(9);
2999 
3000  // f - flag, f=0 for determination of Hessian and gradient of the functional,
3001  // f=1 for determination of the functional value only;
3002  // adaptivity is determined on the first step for adp>0 (nodal based)
3003  if (f == 0)
3004  {
3005  if (adp > 0)
3006  avertex(afun, Gloc, R, cell_in, nvert, adp);
3007  if (adp == 0)
3008  {
3009  for (unsigned i=0; i<_dim; i++)
3010  Gloc[i] = 1.;
3011  }
3012  }
3013 
3014  double
3015  sigma = 0.,
3016  fun = 0,
3017  gqmin = 1e32,
3018  g = 0, // Vmin
3019  lqmin = 0.;
3020 
3021  // cell integration depending on cell type
3022  if (_dim == 2)
3023  {
3024  // 2D
3025  if (nvert == 3)
3026  {
3027  // tri
3028  sigma = 1.;
3029  fun += vertex(W, F, R, cell_in, epsilon, w, nvert, K, H, me, vol, f, lqmin, adp, Gloc, sigma);
3030  g += sigma*lqmin;
3031  if (gqmin > lqmin)
3032  gqmin = lqmin;
3033  }
3034  if (nvert == 4)
3035  {
3036  //quad
3037  for (unsigned i=0; i<2; i++)
3038  {
3039  K[0] = i;
3040  for (unsigned j=0; j<2; j++)
3041  {
3042  K[1] = j;
3043  sigma = 0.25;
3044  fun += vertex(W, F, R, cell_in, epsilon, w, nvert, K, H, me,
3045  vol, f, lqmin, adp, Gloc, sigma);
3046  g += sigma*lqmin;
3047  if (gqmin > lqmin)
3048  gqmin = lqmin;
3049  }
3050  }
3051  }
3052  else
3053  {
3054  // quad tri
3055  for (unsigned i=0; i<3; i++)
3056  {
3057  K[0] = i*0.5;
3058  unsigned k = i/2;
3059  K[1] = static_cast<double>(k);
3060 
3061  for (unsigned j=0; j<3; j++)
3062  {
3063  K[2] = j*0.5;
3064  k = j/2;
3065  K[3] = static_cast<double>(k);
3066  if (i == j)
3067  sigma = 1./12.;
3068  else
3069  sigma = 1./24.;
3070 
3071  fun += vertex(W, F, R, cell_in, epsilon, w, nvert, K, H, me,
3072  vol, f, lqmin, adp, Gloc, sigma);
3073  g += sigma*lqmin;
3074  if (gqmin > lqmin)
3075  gqmin = lqmin;
3076  }
3077  }
3078  }
3079  }
3080  if (_dim == 3)
3081  {
3082  // 3D
3083  if (nvert == 4)
3084  {
3085  // tetr
3086  sigma = 1.;
3087  fun += vertex(W, F, R, cell_in, epsilon, w, nvert, K, H, me,
3088  vol, f, lqmin, adp, Gloc, sigma);
3089  g += sigma*lqmin;
3090  if (gqmin > lqmin)
3091  gqmin = lqmin;
3092  }
3093  if (nvert == 6)
3094  {
3095  //prism
3096  for (unsigned i=0; i<2; i++)
3097  {
3098  K[0] = i;
3099  for (unsigned j=0; j<2; j++)
3100  {
3101  K[1] = j;
3102  for (unsigned k=0; k<3; k++)
3103  {
3104  K[2] = 0.5*static_cast<double>(k);
3105  K[3] = static_cast<double>(k % 2);
3106  sigma = 1./12.;
3107  fun += vertex(W, F, R, cell_in, epsilon, w, nvert, K, H, me,
3108  vol, f, lqmin, adp, Gloc, sigma);
3109  g += sigma*lqmin;
3110  if (gqmin > lqmin)
3111  gqmin = lqmin;
3112  }
3113  }
3114  }
3115  }
3116  if (nvert == 8)
3117  {
3118  // hex
3119  for (unsigned i=0; i<2; i++)
3120  {
3121  K[0] = i;
3122  for (unsigned j=0; j<2; j++)
3123  {
3124  K[1] = j;
3125  for (unsigned k=0; k<2; k++)
3126  {
3127  K[2] = k;
3128  for (unsigned l=0; l<2; l++)
3129  {
3130  K[3] = l;
3131  for (unsigned m=0; m<2; m++)
3132  {
3133  K[4] = m;
3134  for (unsigned nn=0; nn<2; nn++)
3135  {
3136  K[5] = nn;
3137 
3138  if ((i==nn) && (j==l) && (k==m))
3139  sigma = 1./27.;
3140 
3141  if (((i==nn) && (j==l) && (k!=m)) ||
3142  ((i==nn) && (j!=l) && (k==m)) ||
3143  ((i!=nn) && (j==l) && (k==m)))
3144  sigma = 1./54.;
3145 
3146  if (((i==nn) && (j!=l) && (k!=m)) ||
3147  ((i!=nn) && (j!=l) && (k==m)) ||
3148  ((i!=nn) && (j==l) && (k!=m)))
3149  sigma = 1./108.;
3150 
3151  if ((i!=nn) && (j!=l) && (k!=m))
3152  sigma=1./216.;
3153 
3154  fun += vertex(W, F, R, cell_in, epsilon, w, nvert, K, H, me,
3155  vol, f, lqmin, adp, Gloc, sigma);
3156  g += sigma*lqmin;
3157 
3158  if (gqmin > lqmin)
3159  gqmin = lqmin;
3160  }
3161  }
3162  }
3163  }
3164  }
3165  }
3166  }
3167  else
3168  {
3169  // quad tetr
3170  for (unsigned i=0; i<4; i++)
3171  {
3172  for (unsigned j=0; j<4; j++)
3173  {
3174  for (unsigned k=0; k<4; k++)
3175  {
3176  switch (i)
3177  {
3178  case 0:
3179  K[0] = 0;
3180  K[1] = 0;
3181  K[2] = 0;
3182  break;
3183 
3184  case 1:
3185  K[0] = 1;
3186  K[1] = 0;
3187  K[2] = 0;
3188  break;
3189 
3190  case 2:
3191  K[0] = 0.5;
3192  K[1] = 1;
3193  K[2] = 0;
3194  break;
3195 
3196  case 3:
3197  K[0] = 0.5;
3198  K[1] = 1./3.;
3199  K[2] = 1;
3200  break;
3201 
3202  default:
3203  break;
3204  }
3205 
3206  switch (j)
3207  {
3208  case 0:
3209  K[3] = 0;
3210  K[4] = 0;
3211  K[5] = 0;
3212  break;
3213 
3214  case 1:
3215  K[3] = 1;
3216  K[4] = 0;
3217  K[5] = 0;
3218  break;
3219 
3220  case 2:
3221  K[3] = 0.5;
3222  K[4] = 1;
3223  K[5] = 0;
3224  break;
3225 
3226  case 3:
3227  K[3] = 0.5;
3228  K[4] = 1./3.;
3229  K[5] = 1;
3230  break;
3231 
3232  default:
3233  break;
3234  }
3235 
3236  switch (k)
3237  {
3238  case 0:
3239  K[6] = 0;
3240  K[7] = 0;
3241  K[8] = 0;
3242  break;
3243 
3244  case 1:
3245  K[6] = 1;
3246  K[7] = 0;
3247  K[8] = 0;
3248  break;
3249 
3250  case 2:
3251  K[6] = 0.5;
3252  K[7] = 1;
3253  K[8] = 0;
3254  break;
3255 
3256  case 3:
3257  K[6] = 0.5;
3258  K[7] = 1./3.;
3259  K[8] = 1;
3260  break;
3261 
3262  default:
3263  break;
3264  }
3265 
3266  if ((i==j) && (j==k))
3267  sigma = 1./120.;
3268 
3269  else if ((i==j) || (j==k) || (i==k))
3270  sigma = 1./360.;
3271 
3272  else
3273  sigma = 1./720.;
3274 
3275  fun += vertex(W, F, R, cell_in, epsilon, w, nvert, K, H, me,
3276  vol, f, lqmin, adp, Gloc, sigma);
3277 
3278  g += sigma*lqmin;
3279  if (gqmin > lqmin)
3280  gqmin = lqmin;
3281  }
3282  }
3283  }
3284  }
3285  }
3286 
3287  // fixed nodes correction
3288  for (int ii=0; ii<nvert; ii++)
3289  {
3290  if (mask[cell_in[ii]] == 1)
3291  {
3292  for (unsigned kk=0; kk<_dim; kk++)
3293  {
3294  for (int jj=0; jj<nvert; jj++)
3295  {
3296  W[kk][ii][jj] = 0;
3297  W[kk][jj][ii] = 0;
3298  }
3299 
3300  W[kk][ii][ii] = 1;
3301  F[kk][ii] = 0;
3302  }
3303  }
3304  }
3305  // end of fixed nodes correction
3306 
3307  // Set up return value references
3308  Vmin = g;
3309  qmin = gqmin/vol;
3310 
3311  return fun;
3312 }
3313 
3314 
3315 
3316 // avertex - assembly of adaptivity metric on a cell
3317 double VariationalMeshSmoother::avertex(const std::vector<double> & afun,
3318  std::vector<double> & G,
3319  const Array2D<double> & R,
3320  const std::vector<int> & cell_in,
3321  int nvert,
3322  int adp)
3323 {
3324  std::vector<double> a1(3), a2(3), a3(3), qu(3), K(8);
3325  Array2D<double> Q(3, nvert);
3326 
3327  for (int i=0; i<8; i++)
3328  K[i] = 0.5; // cell center
3329 
3330  basisA(Q, nvert, K, Q, 1);
3331 
3332  for (unsigned i=0; i<_dim; i++)
3333  for (int j=0; j<nvert; j++)
3334  {
3335  a1[i] += Q[i][j]*R[cell_in[j]][0];
3336  a2[i] += Q[i][j]*R[cell_in[j]][1];
3337  if (_dim == 3)
3338  a3[i] += Q[i][j]*R[cell_in[j]][2];
3339 
3340  qu[i] += Q[i][j]*afun[cell_in[j]];
3341  }
3342 
3343  double det = 0.;
3344 
3345  if (_dim == 2)
3346  det = jac2(a1[0], a1[1],
3347  a2[0], a2[1]);
3348  else
3349  det = jac3(a1[0], a1[1], a1[2],
3350  a2[0], a2[1], a2[2],
3351  a3[0], a3[1], a3[2]);
3352 
3353  // Return val
3354  double g = 0.;
3355 
3356  if (det != 0)
3357  {
3358  if (_dim == 2)
3359  {
3360  double df0 = jac2(qu[0], qu[1], a2[0], a2[1])/det;
3361  double df1 = jac2(a1[0], a1[1], qu[0], qu[1])/det;
3362  g = (1. + df0*df0 + df1*df1);
3363 
3364  if (adp == 2)
3365  {
3366  // directional adaptation G=diag(g_i)
3367  G[0] = 1. + df0*df0;
3368  G[1] = 1. + df1*df1;
3369  }
3370  else
3371  {
3372  for (unsigned i=0; i<_dim; i++)
3373  G[i] = g; // simple adaptation G=gI
3374  }
3375  }
3376  else
3377  {
3378  double df0 = (qu[0]*(a2[1]*a3[2]-a2[2]*a3[1]) + qu[1]*(a2[2]*a3[0]-a2[0]*a3[2]) + qu[2]*(a2[0]*a3[1]-a2[1]*a3[0]))/det;
3379  double df1 = (qu[0]*(a3[1]*a1[2]-a3[2]*a1[1]) + qu[1]*(a3[2]*a1[0]-a3[0]*a1[2]) + qu[2]*(a3[0]*a1[1]-a3[1]*a1[0]))/det;
3380  double df2 = (qu[0]*(a1[1]*a2[2]-a1[2]*a2[1]) + qu[1]*(a1[2]*a2[0]-a1[0]*a2[2]) + qu[2]*(a1[0]*a2[1]-a1[1]*a2[0]))/det;
3381  g = 1. + df0*df0 + df1*df1 + df2*df2;
3382  if (adp == 2)
3383  {
3384  // directional adaptation G=diag(g_i)
3385  G[0] = 1. + df0*df0;
3386  G[1] = 1. + df1*df1;
3387  G[2] = 1. + df2*df2;
3388  }
3389  else
3390  {
3391  for (unsigned i=0; i<_dim; i++)
3392  G[i] = g; // simple adaptation G=gI
3393  }
3394  }
3395  }
3396  else
3397  {
3398  g = 1.;
3399  for (unsigned i=0; i<_dim; i++)
3400  G[i] = g;
3401  }
3402 
3403  return g;
3404 }
3405 
3406 
3407 
3408 // Computes local matrix W and local rhs F on one basis
3409 double VariationalMeshSmoother::vertex(Array3D<double> & W,
3410  Array2D<double> & F,
3411  const Array2D<double> & R,
3412  const std::vector<int> & cell_in,
3413  double epsilon,
3414  double w,
3415  int nvert,
3416  const std::vector<double> & K,
3417  const Array2D<double> & H,
3418  int me,
3419  double vol,
3420  int f,
3421  double & qmin,
3422  int adp,
3423  const std::vector<double> & g,
3424  double sigma)
3425 {
3426  // hessian, function, gradient
3427  Array2D<double> Q(3, nvert);
3428  basisA(Q, nvert, K, H, me);
3429 
3430  std::vector<double> a1(3), a2(3), a3(3);
3431  for (unsigned i=0; i<_dim; i++)
3432  for (int j=0; j<nvert; j++)
3433  {
3434  a1[i] += Q[i][j]*R[cell_in[j]][0];
3435  a2[i] += Q[i][j]*R[cell_in[j]][1];
3436  if (_dim == 3)
3437  a3[i] += Q[i][j]*R[cell_in[j]][2];
3438  }
3439 
3440  // account for adaptation
3441  double G = 0.;
3442  if (adp < 2)
3443  G = g[0];
3444  else
3445  {
3446  G = 1.;
3447  for (unsigned i=0; i<_dim; i++)
3448  {
3449  a1[i] *= std::sqrt(g[0]);
3450  a2[i] *= std::sqrt(g[1]);
3451  if (_dim == 3)
3452  a3[i] *= std::sqrt(g[2]);
3453  }
3454  }
3455 
3456  double
3457  det = 0.,
3458  tr = 0.,
3459  phit = 0.;
3460 
3461  std::vector<double> av1(3), av2(3), av3(3);
3462 
3463  if (_dim == 2)
3464  {
3465  av1[0] = a2[1];
3466  av1[1] = -a2[0];
3467  av2[0] = -a1[1];
3468  av2[1] = a1[0];
3469  det = jac2(a1[0], a1[1], a2[0], a2[1]);
3470  for (int i=0; i<2; i++)
3471  tr += 0.5*(a1[i]*a1[i] + a2[i]*a2[i]);
3472 
3473  phit = (1-w)*tr + w*0.5*(vol/G + det*det*G/vol);
3474  }
3475 
3476  if (_dim == 3)
3477  {
3478  av1[0] = a2[1]*a3[2] - a2[2]*a3[1];
3479  av1[1] = a2[2]*a3[0] - a2[0]*a3[2];
3480  av1[2] = a2[0]*a3[1] - a2[1]*a3[0];
3481 
3482  av2[0] = a3[1]*a1[2] - a3[2]*a1[1];
3483  av2[1] = a3[2]*a1[0] - a3[0]*a1[2];
3484  av2[2] = a3[0]*a1[1] - a3[1]*a1[0];
3485 
3486  av3[0] = a1[1]*a2[2] - a1[2]*a2[1];
3487  av3[1] = a1[2]*a2[0] - a1[0]*a2[2];
3488  av3[2] = a1[0]*a2[1] - a1[1]*a2[0];
3489 
3490  det = jac3(a1[0], a1[1], a1[2],
3491  a2[0], a2[1], a2[2],
3492  a3[0], a3[1], a3[2]);
3493 
3494  for (int i=0; i<3; i++)
3495  tr += (1./3.)*(a1[i]*a1[i] + a2[i]*a2[i] + a3[i]*a3[i]);
3496 
3497  phit = (1-w)*pow(tr,1.5) + w*0.5*(vol/G + det*det*G/vol);
3498  }
3499 
3500  double dchi = 0.5 + det*0.5/std::sqrt(epsilon*epsilon + det*det);
3501  double chi = 0.5*(det + std::sqrt(epsilon*epsilon + det*det));
3502  double fet = (chi != 0.) ? phit/chi : 1.e32;
3503 
3504  // Set return value reference
3505  qmin = det*G;
3506 
3507  // gradient and Hessian
3508  if (f == 0)
3509  {
3510  Array3D<double> P(3, 3, 3);
3511  Array3D<double> d2phi(3, 3, 3);
3512  Array2D<double> dphi(3, 3);
3513  Array2D<double> dfe(3, 3);
3514 
3515  if (_dim == 2)
3516  {
3517  for (int i=0; i<2; i++)
3518  {
3519  dphi[0][i] = (1-w)*a1[i] + w*G*det*av1[i]/vol;
3520  dphi[1][i] = (1-w)*a2[i] + w*G*det*av2[i]/vol;
3521  }
3522 
3523  for (int i=0; i<2; i++)
3524  for (int j=0; j<2; j++)
3525  {
3526  d2phi[0][i][j] = w*G*av1[i]*av1[j]/vol;
3527  d2phi[1][i][j] = w*G*av2[i]*av2[j]/vol;
3528 
3529  if (i == j)
3530  for (int k=0; k<2; k++)
3531  d2phi[k][i][j] += 1.-w;
3532  }
3533 
3534  for (int i=0; i<2; i++)
3535  {
3536  dfe[0][i] = (dphi[0][i] - fet*dchi*av1[i])/chi;
3537  dfe[1][i] = (dphi[1][i] - fet*dchi*av2[i])/chi;
3538  }
3539 
3540  for (int i=0; i<2; i++)
3541  for (int j=0; j<2; j++)
3542  {
3543  P[0][i][j] = (d2phi[0][i][j] - dchi*(av1[i]*dfe[0][j] + dfe[0][i]*av1[j]))/chi;
3544  P[1][i][j] = (d2phi[1][i][j] - dchi*(av2[i]*dfe[1][j] + dfe[1][i]*av2[j]))/chi;
3545  }
3546  }
3547 
3548  if (_dim == 3)
3549  {
3550  for (int i=0; i<3; i++)
3551  {
3552  dphi[0][i] = (1-w)*std::sqrt(tr)*a1[i] + w*G*det*av1[i]/vol;
3553  dphi[1][i] = (1-w)*std::sqrt(tr)*a2[i] + w*G*det*av2[i]/vol;
3554  dphi[2][i] = (1-w)*std::sqrt(tr)*a3[i] + w*G*det*av3[i]/vol;
3555  }
3556 
3557  for (int i=0; i<3; i++)
3558  {
3559  if (tr != 0)
3560  {
3561  d2phi[0][i][i] = (1-w)/(3.*std::sqrt(tr))*a1[i]*a1[i] + w*G*av1[i]*av1[i]/vol;
3562  d2phi[1][i][i] = (1-w)/(3.*std::sqrt(tr))*a2[i]*a2[i] + w*G*av2[i]*av2[i]/vol;
3563  d2phi[2][i][i] = (1-w)/(3.*std::sqrt(tr))*a3[i]*a3[i] + w*G*av3[i]*av3[i]/vol;
3564  }
3565  else
3566  {
3567  for (int k=0; k<3; k++)
3568  d2phi[k][i][i] = 0.;
3569  }
3570 
3571  for (int k=0; k<3; k++)
3572  d2phi[k][i][i] += (1-w)*std::sqrt(tr);
3573  }
3574 
3575  const double con = 100.;
3576 
3577  for (int i=0; i<3; i++)
3578  {
3579  dfe[0][i] = (dphi[0][i] - fet*dchi*av1[i] + con*epsilon*a1[i])/chi;
3580  dfe[1][i] = (dphi[1][i] - fet*dchi*av2[i] + con*epsilon*a2[i])/chi;
3581  dfe[2][i] = (dphi[2][i] - fet*dchi*av3[i] + con*epsilon*a3[i])/chi;
3582  }
3583 
3584  for (int i=0; i<3; i++)
3585  {
3586  P[0][i][i] = (d2phi[0][i][i] - dchi*(av1[i]*dfe[0][i] + dfe[0][i]*av1[i]) + con*epsilon)/chi;
3587  P[1][i][i] = (d2phi[1][i][i] - dchi*(av2[i]*dfe[1][i] + dfe[1][i]*av2[i]) + con*epsilon)/chi;
3588  P[2][i][i] = (d2phi[2][i][i] - dchi*(av3[i]*dfe[2][i] + dfe[2][i]*av3[i]) + con*epsilon)/chi;
3589  }
3590 
3591  for (int k=0; k<3; k++)
3592  for (int i=0; i<3; i++)
3593  for (int j=0; j<3; j++)
3594  if (i != j)
3595  P[k][i][j] = 0.;
3596  }
3597 
3598  /*--------matrix W, right side F---------------------*/
3599  Array2D<double> gpr(3, nvert);
3600 
3601  for (unsigned i1=0; i1<_dim; i1++)
3602  {
3603  for (unsigned i=0; i<_dim; i++)
3604  for (int j=0; j<nvert; j++)
3605  for (unsigned k=0; k<_dim; k++)
3606  gpr[i][j] += P[i1][i][k]*Q[k][j];
3607 
3608  for (int i=0; i<nvert; i++)
3609  for (int j=0; j<nvert; j++)
3610  for (unsigned k=0; k<_dim; k++)
3611  W[i1][i][j] += Q[k][i]*gpr[k][j]*sigma;
3612 
3613  for (int i=0; i<nvert; i++)
3614  for (int k=0; k<2; k++)
3615  F[i1][i] += Q[k][i]*dfe[i1][k]*sigma;
3616  }
3617  }
3618 
3619  return fet*sigma;
3620 }
3621 
3622 
3623 
3624 // Solve Symmetric Positive Definite (SPD) linear system
3625 // by Conjugate Gradient (CG) method preconditioned by
3626 // Point Jacobi (diagonal scaling)
3627 int VariationalMeshSmoother::solver(int n,
3628  const std::vector<int> & ia,
3629  const std::vector<int> & ja,
3630  const std::vector<double> & a,
3631  std::vector<double> & x,
3632  const std::vector<double> & b,
3633  double eps,
3634  int maxite,
3635  int msglev)
3636 {
3637  _logfile << "Beginning Solve " << n << std::endl;
3638 
3639  // Check parameters
3640  int info = pcg_par_check(n, ia, ja, a, eps, maxite, msglev);
3641  if (info != 0)
3642  return info;
3643 
3644  // PJ preconditioner construction
3645  std::vector<double> u(n);
3646  for (int i=0; i<n; i++)
3647  u[i] = 1./a[ia[i]];
3648 
3649  // PCG iterations
3650  std::vector<double> r(n), p(n), z(n);
3651  info = pcg_ic0(n, ia, ja, a, u, x, b, r, p, z, eps, maxite, msglev);
3652 
3653  // Mat sparse_global;
3654  // MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,n,n,ia,ja,a,&sparse_global);
3655 
3656  _logfile << "Finished Solve " << n << std::endl;
3657 
3658  return info;
3659 }
3660 
3661 
3662 
3663 // Input parameter check
3664 int VariationalMeshSmoother::pcg_par_check(int n,
3665  const std::vector<int> & ia,
3666  const std::vector<int> & ja,
3667  const std::vector<double> & a,
3668  double eps,
3669  int maxite,
3670  int msglev)
3671 {
3672  int i, j, jj, k;
3673 
3674  if (n <= 0)
3675  {
3676  if (msglev > 0)
3677  _logfile << "sol_pcg: incorrect input parameter: n =" << n << std::endl;
3678  return -1;
3679  }
3680 
3681  if (ia[0] != 0)
3682  {
3683  if (msglev > 0)
3684  _logfile << "sol_pcg: incorrect input parameter: ia(1) =" << ia[0] << std::endl;
3685  return -2;
3686  }
3687 
3688  for (i=0; i<n; i++)
3689  {
3690  if (ia[i+1] < ia[i])
3691  {
3692  if (msglev >= 1)
3693  _logfile << "sol_pcg: incorrect input parameter: i ="
3694  << i+1
3695  << "; ia(i) ="
3696  << ia[i]
3697  << " must be <= a(i+1) ="
3698  << ia[i+1]
3699  << std::endl;
3700  return -2;
3701  }
3702  }
3703 
3704  for (i=0; i<n; i++)
3705  {
3706  if (ja[ia[i]] != (i+1))
3707  {
3708  if (msglev >= 1)
3709  _logfile << "sol_pcg: incorrect input parameter: i ="
3710  << i+1
3711  << " ; ia(i) ="
3712  << ia[i]
3713  << " ; absence of diagonal entry; k ="
3714  << ia[i]+1
3715  << " ; the value ja(k) ="
3716  << ja[ia[i]]
3717  << " must be equal to i"
3718  << std::endl;
3719 
3720  return -3;
3721  }
3722 
3723  jj = 0;
3724  for (k=ia[i]; k<ia[i+1]; k++)
3725  {
3726  j=ja[k];
3727  if ((j<(i+1)) || (j>n))
3728  {
3729  if (msglev >= 1)
3730  _logfile << "sol_pcg: incorrect input parameter: i ="
3731  << i+1
3732  << " ; ia(i) ="
3733  << ia[i]
3734  << " ; a(i+1) ="
3735  << ia[i+1]
3736  << " ; k ="
3737  << k+1
3738  << " ; the value ja(k) ="
3739  << ja[k]
3740  << " is out of range"
3741  << std::endl;
3742 
3743  return -3;
3744  }
3745  if (j <= jj)
3746  {
3747  if (msglev >= 1)
3748  _logfile << "sol_pcg: incorrect input parameter: i ="
3749  << i+1
3750  << " ; ia(i) ="
3751  << ia[i]
3752  << " ; a(i+1) ="
3753  << ia[i+1]
3754  << " ; k ="
3755  << k+1
3756  << " ; the value ja(k) ="
3757  << ja[k]
3758  << " must be less than ja(k+1) ="
3759  << ja[k+1]
3760  << std::endl;
3761 
3762  return -3;
3763  }
3764  jj = j;
3765  }
3766  }
3767 
3768  for (i=0; i<n; i++)
3769  {
3770  if (a[ia[i]] <= 0.)
3771  {
3772  if (msglev >= 1)
3773  _logfile << "sol_pcg: incorrect input parameter: i ="
3774  << i+1
3775  << " ; ia(i) ="
3776  << ia[i]
3777  << " ; the diagonal entry a(ia(i)) ="
3778  << a[ia[i]]
3779  << " must be > 0"
3780  << std::endl;
3781 
3782  return -4;
3783  }
3784  }
3785 
3786  if (eps < 0.)
3787  {
3788  if (msglev > 0)
3789  _logfile << "sol_pcg: incorrect input parameter: eps =" << eps << std::endl;
3790  return -7;
3791  }
3792 
3793  if (maxite < 1)
3794  {
3795  if (msglev > 0)
3796  _logfile << "sol_pcg: incorrect input parameter: maxite =" << maxite << std::endl;
3797  return -8;
3798  }
3799 
3800  return 0;
3801 }
3802 
3803 
3804 
3805 // Solve the SPD linear system by PCG method
3806 int VariationalMeshSmoother::pcg_ic0(int n,
3807  const std::vector<int> & ia,
3808  const std::vector<int> & ja,
3809  const std::vector<double> & a,
3810  const std::vector<double> & u,
3811  std::vector<double> & x,
3812  const std::vector<double> & b,
3813  std::vector<double> & r,
3814  std::vector<double> & p,
3815  std::vector<double> & z,
3816  double eps,
3817  int maxite,
3818  int msglev)
3819 {
3820  // Return value
3821  int i = 0;
3822 
3823  double
3824  rhr = 0.,
3825  rhr0 = 0.,
3826  rhrold = 0.,
3827  rr0 = 0.;
3828 
3829  for (i=0; i<=maxite; i++)
3830  {
3831  if (i == 0)
3832  p = x;
3833 
3834  // z=Ap
3835  for (int ii=0; ii<n; ii++)
3836  z[ii] = 0.;
3837 
3838  for (int ii=0; ii<n; ii++)
3839  {
3840  z[ii] += a[ia[ii]]*p[ii];
3841 
3842  for (int j=ia[ii]+1; j<ia[ii+1]; j++)
3843  {
3844  z[ii] += a[j]*p[ja[j]-1];
3845  z[ja[j]-1] += a[j]*p[ii];
3846  }
3847  }
3848 
3849  if (i == 0)
3850  for (int k=0; k<n; k++)
3851  r[k] = b[k] - z[k]; // r(0) = b - Ax(0)
3852 
3853  if (i > 0)
3854  {
3855  double pap = 0.;
3856  for (int k=0; k<n; k++)
3857  pap += p[k]*z[k];
3858 
3859  double alpha = rhr/pap;
3860  for (int k=0; k<n; k++)
3861  {
3862  x[k] += p[k]*alpha;
3863  r[k] -= z[k]*alpha;
3864  }
3865  rhrold = rhr;
3866  }
3867 
3868  // z = H * r
3869  for (int ii=0; ii<n; ii++)
3870  z[ii] = r[ii]*u[ii];
3871 
3872  if (i == 0)
3873  p = z;
3874 
3875  rhr = 0.;
3876  double rr = 0.;
3877  for (int k=0; k<n; k++)
3878  {
3879  rhr += r[k]*z[k];
3880  rr += r[k]*r[k];
3881  }
3882 
3883  if (i == 0)
3884  {
3885  rhr0 = rhr;
3886  rr0 = rr;
3887  }
3888 
3889  if (msglev > 0)
3890  _logfile << i
3891  << " ) rHr ="
3892  << std::sqrt(rhr/rhr0)
3893  << " rr ="
3894  << std::sqrt(rr/rr0)
3895  << std::endl;
3896 
3897  if (rr <= eps*eps*rr0)
3898  return i;
3899 
3900  if (i >= maxite)
3901  return i;
3902 
3903  if (i > 0)
3904  {
3905  double beta = rhr/rhrold;
3906  for (int k=0; k<n; k++)
3907  p[k] = z[k] + p[k]*beta;
3908  }
3909  } // end for i
3910 
3911  return i;
3912 }
3913 
3914 
3915 
3916 
3917 // Sample mesh file generation
3918 void VariationalMeshSmoother::gener(char grid[],
3919  int n)
3920 {
3921  const int n1 = 3;
3922 
3923  int
3924  N = 1,
3925  ncells = 1;
3926 
3927  for (int i=0; i<n; i++)
3928  {
3929  N *= n1;
3930  ncells *= (n1-1);
3931  }
3932 
3933  const double x = 1./static_cast<double>(n1-1);
3934 
3935  std::ofstream outfile(grid);
3936 
3937  outfile << n << "\n" << N << "\n" << ncells << "\n0\n" << std::endl;
3938 
3939  for (int i=0; i<N; i++)
3940  {
3941  // node coordinates
3942  int k = i;
3943  int mask = 0;
3944  for (int j=0; j<n; j++)
3945  {
3946  int ii = k % n1;
3947  if ((ii == 0) || (ii == n1-1))
3948  mask = 1;
3949 
3950  outfile << static_cast<double>(ii)*x << " ";
3951  k /= n1;
3952  }
3953  outfile << mask << std::endl;
3954  }
3955 
3956  for (int i=0; i<ncells; i++)
3957  {
3958  // cell connectivity
3959  int nc = i;
3960  int ii = nc%(n1-1);
3961  nc /= (n1-1);
3962  int jj = nc%(n1-1);
3963  int kk = nc/(n1-1);
3964 
3965  if (n == 2)
3966  outfile << ii+n1*jj << " "
3967  << ii+1+jj*n1 << " "
3968  << ii+(jj+1)*n1 << " "
3969  << ii+1+(jj+1)*n1 << " ";
3970 
3971  if (n == 3)
3972  outfile << ii+n1*jj+n1*n1*kk << " "
3973  << ii+1+jj*n1+n1*n1*kk << " "
3974  << ii+(jj+1)*n1+n1*n1*kk << " "
3975  << ii+1+(jj+1)*n1+n1*n1*kk << " "
3976  << ii+n1*jj+n1*n1*(kk+1) << " "
3977  << ii+1+jj*n1+n1*n1*(kk+1) << " "
3978  << ii+(jj+1)*n1+n1*n1*(kk+1) << " "
3979  << ii+1+(jj+1)*n1+n1*n1*(kk+1) << " ";
3980 
3981  outfile << "-1 -1 0 \n";
3982  }
3983 }
3984 
3985 
3986 
3987 // Metric Generation
3988 void VariationalMeshSmoother::metr_data_gen(std::string grid,
3989  std::string metr,
3990  int me)
3991 {
3992  double det, g1, g2, g3, det_o, g1_o, g2_o, g3_o, eps=1e-3;
3993 
3994  std::vector<double> K(9);
3995  Array2D<double> Q(3, 3*_dim + _dim%2);
3996 
3997  // Use _mesh to determine N and ncells
3998  this->_n_nodes = _mesh.n_nodes();
3999  this->_n_cells = _mesh.n_active_elem();
4000 
4001  std::vector<int>
4002  mask(_n_nodes),
4003  mcells(_n_cells);
4004 
4005  Array2D<int> cells(_n_cells, 3*_dim + _dim%2);
4006  Array2D<double> R(_n_nodes,_dim);
4007 
4008  readgr(R, mask, cells, mcells, mcells, mcells);
4009 
4010  // generate metric file
4011  std::ofstream metric_file(metr.c_str());
4012  if (!metric_file.good())
4013  libmesh_error_msg("Error opening metric output file.");
4014 
4015  // Use scientific notation with 6 digits
4016  metric_file << std::scientific << std::setprecision(6);
4017 
4018  int Ncells = 0;
4019  det_o = 1.;
4020  g1_o = 1.;
4021  g2_o = 1.;
4022  g3_o = 1.;
4023  for (unsigned i=0; i<_n_cells; i++)
4024  if (mcells[i] >= 0)
4025  {
4026  int nvert=0;
4027  while (cells[i][nvert] >= 0)
4028  nvert++;
4029 
4030  if (_dim == 2)
4031  {
4032  // 2D - tri and quad
4033  if (nvert == 3)
4034  {
4035  // tri
4036  basisA(Q, 3, K, Q, 1);
4037 
4038  Point a1, a2;
4039  for (int k=0; k<2; k++)
4040  for (int l=0; l<3; l++)
4041  {
4042  a1(k) += Q[k][l]*R[cells[i][l]][0];
4043  a2(k) += Q[k][l]*R[cells[i][l]][1];
4044  }
4045 
4046  det = jac2(a1(0), a1(1), a2(0), a2(1));
4047  g1 = std::sqrt(a1(0)*a1(0) + a2(0)*a2(0));
4048  g2 = std::sqrt(a1(1)*a1(1) + a2(1)*a2(1));
4049 
4050  // need to keep data from previous cell
4051  if ((std::abs(det) < eps*eps*det_o) || (det < 0))
4052  det = det_o;
4053 
4054  if ((std::abs(g1) < eps*g1_o) || (g1<0))
4055  g1 = g1_o;
4056 
4057  if ((std::abs(g2) < eps*g2_o) || (g2<0))
4058  g2 = g2_o;
4059 
4060  // write to file
4061  if (me == 2)
4062  metric_file << 1./std::sqrt(det)
4063  << " 0.000000e+00 \n0.000000e+00 "
4064  << 1./std::sqrt(det)
4065  << std::endl;
4066 
4067  if (me == 3)
4068  metric_file << 1./g1
4069  << " 0.000000e+00 \n0.000000e+00 "
4070  << 1./g2
4071  << std::endl;
4072 
4073  det_o = det;
4074  g1_o = g1;
4075  g2_o = g2;
4076  Ncells++;
4077  }
4078 
4079  if (nvert == 4)
4080  {
4081  // quad
4082 
4083  // Set up the node edge neighbor indices
4084  const unsigned node_indices[4] = {0, 1, 3, 2};
4085  const unsigned first_neighbor_indices[4] = {1, 3, 2, 0};
4086  const unsigned second_neighbor_indices[4] = {2, 0, 1, 3};
4087 
4088  // Loop over each node, compute some quantities associated
4089  // with its edge neighbors, and write them to file.
4090  for (unsigned ni=0; ni<4; ++ni)
4091  {
4092  unsigned
4093  node_index = node_indices[ni],
4094  first_neighbor_index = first_neighbor_indices[ni],
4095  second_neighbor_index = second_neighbor_indices[ni];
4096 
4097  Real
4098  node_x = R[cells[i][node_index]][0],
4099  node_y = R[cells[i][node_index]][1],
4100  first_neighbor_x = R[cells[i][first_neighbor_index]][0],
4101  first_neighbor_y = R[cells[i][first_neighbor_index]][1],
4102  second_neighbor_x = R[cells[i][second_neighbor_index]][0],
4103  second_neighbor_y = R[cells[i][second_neighbor_index]][1];
4104 
4105 
4106  // "dx"
4107  Point a1(first_neighbor_x - node_x,
4108  second_neighbor_x - node_x);
4109 
4110  // "dy"
4111  Point a2(first_neighbor_y - node_y,
4112  second_neighbor_y - node_y);
4113 
4114  det = jac2(a1(0), a1(1), a2(0), a2(1));
4115  g1 = std::sqrt(a1(0)*a1(0) + a2(0)*a2(0));
4116  g2 = std::sqrt(a1(1)*a1(1) + a2(1)*a2(1));
4117 
4118  // need to keep data from previous cell
4119  if ((std::abs(det) < eps*eps*det_o) || (det < 0))
4120  det = det_o;
4121 
4122  if ((std::abs(g1) < eps*g1_o) || (g1 < 0))
4123  g1 = g1_o;
4124 
4125  if ((std::abs(g2) < eps*g2_o) || (g2 < 0))
4126  g2 = g2_o;
4127 
4128  // write to file
4129  if (me == 2)
4130  metric_file << 1./std::sqrt(det) << " "
4131  << 0.5/std::sqrt(det) << " \n0.000000e+00 "
4132  << 0.5*std::sqrt(3./det)
4133  << std::endl;
4134 
4135  if (me == 3)
4136  metric_file << 1./g1 << " "
4137  << 0.5/g2 << "\n0.000000e+00 "
4138  << 0.5*std::sqrt(3.)/g2
4139  << std::endl;
4140 
4141  det_o = det;
4142  g1_o = g1;
4143  g2_o = g2;
4144  Ncells++;
4145  }
4146  } // end QUAD case
4147  } // end _dim==2
4148 
4149  if (_dim == 3)
4150  {
4151  // 3D - tetr and hex
4152 
4153  if (nvert == 4)
4154  {
4155  // tetr
4156  basisA(Q, 4, K, Q, 1);
4157 
4158  Point a1, a2, a3;
4159 
4160  for (int k=0; k<3; k++)
4161  for (int l=0; l<4; l++)
4162  {
4163  a1(k) += Q[k][l]*R[cells[i][l]][0];
4164  a2(k) += Q[k][l]*R[cells[i][l]][1];
4165  a3(k) += Q[k][l]*R[cells[i][l]][2];
4166  }
4167 
4168  det = jac3(a1(0), a1(1), a1(2),
4169  a2(0), a2(1), a2(2),
4170  a3(0), a3(1), a3(2));
4171  g1 = std::sqrt(a1(0)*a1(0) + a2(0)*a2(0) + a3(0)*a3(0));
4172  g2 = std::sqrt(a1(1)*a1(1) + a2(1)*a2(1) + a3(1)*a3(1));
4173  g3 = std::sqrt(a1(2)*a1(2) + a2(2)*a2(2) + a3(2)*a3(2));
4174 
4175  // need to keep data from previous cell
4176  if ((std::abs(det) < eps*eps*eps*det_o) || (det < 0))
4177  det = det_o;
4178 
4179  if ((std::abs(g1) < eps*g1_o) || (g1 < 0))
4180  g1 = g1_o;
4181 
4182  if ((std::abs(g2) < eps*g2_o) || (g2 < 0))
4183  g2 = g2_o;
4184 
4185  if ((std::abs(g3) < eps*g3_o) || (g3 < 0))
4186  g3 = g3_o;
4187 
4188  // write to file
4189  if (me == 2)
4190  metric_file << 1./pow(det, 1./3.)
4191  << " 0.000000e+00 0.000000e+00\n0.000000e+00 "
4192  << 1./pow(det, 1./3.)
4193  << " 0.000000e+00\n0.000000e+00 0.000000e+00 "
4194  << 1./pow(det, 1./3.)
4195  << std::endl;
4196 
4197  if (me == 3)
4198  metric_file << 1./g1
4199  << " 0.000000e+00 0.000000e+00\n0.000000e+00 "
4200  << 1./g2
4201  << " 0.000000e+00\n0.000000e+00 0.000000e+00 "
4202  << 1./g3
4203  << std::endl;
4204 
4205  det_o = det;
4206  g1_o = g1;
4207  g2_o = g2;
4208  g3_o = g3;
4209  Ncells++;
4210  }
4211 
4212  if (nvert == 8)
4213  {
4214  // hex
4215 
4216  // Set up the node edge neighbor indices
4217  const unsigned node_indices[8] = {0, 1, 3, 2, 4, 5, 7, 6};
4218  const unsigned first_neighbor_indices[8] = {1, 3, 2, 0, 6, 4, 5, 7};
4219  const unsigned second_neighbor_indices[8] = {2, 0, 1, 3, 5, 7, 6, 4};
4220  const unsigned third_neighbor_indices[8] = {4, 5, 7, 6, 0, 1, 3, 2};
4221 
4222  // Loop over each node, compute some quantities associated
4223  // with its edge neighbors, and write them to file.
4224  for (unsigned ni=0; ni<8; ++ni)
4225  {
4226  unsigned
4227  node_index = node_indices[ni],
4228  first_neighbor_index = first_neighbor_indices[ni],
4229  second_neighbor_index = second_neighbor_indices[ni],
4230  third_neighbor_index = third_neighbor_indices[ni];
4231 
4232  Real
4233  node_x = R[cells[i][node_index]][0],
4234  node_y = R[cells[i][node_index]][1],
4235  node_z = R[cells[i][node_index]][2],
4236  first_neighbor_x = R[cells[i][first_neighbor_index]][0],
4237  first_neighbor_y = R[cells[i][first_neighbor_index]][1],
4238  first_neighbor_z = R[cells[i][first_neighbor_index]][2],
4239  second_neighbor_x = R[cells[i][second_neighbor_index]][0],
4240  second_neighbor_y = R[cells[i][second_neighbor_index]][1],
4241  second_neighbor_z = R[cells[i][second_neighbor_index]][2],
4242  third_neighbor_x = R[cells[i][third_neighbor_index]][0],
4243  third_neighbor_y = R[cells[i][third_neighbor_index]][1],
4244  third_neighbor_z = R[cells[i][third_neighbor_index]][2];
4245 
4246  // "dx"
4247  Point a1(first_neighbor_x - node_x,
4248  second_neighbor_x - node_x,
4249  third_neighbor_x - node_x);
4250 
4251  // "dy"
4252  Point a2(first_neighbor_y - node_y,
4253  second_neighbor_y - node_y,
4254  third_neighbor_y - node_y);
4255 
4256  // "dz"
4257  Point a3(first_neighbor_z - node_z,
4258  second_neighbor_z - node_z,
4259  third_neighbor_z - node_z);
4260 
4261  det = jac3(a1(0), a1(1), a1(2),
4262  a2(0), a2(1), a2(2),
4263  a3(0), a3(1), a3(2));
4264  g1 = std::sqrt(a1(0)*a1(0) + a2(0)*a2(0) + a3(0)*a3(0));
4265  g2 = std::sqrt(a1(1)*a1(1) + a2(1)*a2(1) + a3(1)*a3(1));
4266  g3 = std::sqrt(a1(2)*a1(2) + a2(2)*a2(2) + a3(2)*a3(2));
4267 
4268  // need to keep data from previous cell
4269  if ((std::abs(det) < eps*eps*eps*det_o) || (det < 0))
4270  det = det_o;
4271 
4272  if ((std::abs(g1) < eps*g1_o) || (g1 < 0))
4273  g1 = g1_o;
4274 
4275  if ((std::abs(g2) < eps*g2_o) || (g2 < 0))
4276  g2 = g2_o;
4277 
4278  if ((std::abs(g3) < eps*g3_o) || (g3 < 0))
4279  g3 = g3_o;
4280 
4281  // write to file
4282  if (me == 2)
4283  metric_file << 1./pow(det, 1./3.) << " "
4284  << 0.5/pow(det, 1./3.) << " "
4285  << 0.5/pow(det, 1./3.) << "\n0.000000e+00 "
4286  << 0.5*std::sqrt(3.)/pow(det, 1./3.) << " "
4287  << 0.5/(std::sqrt(3.)*pow(det, 1./3.)) << "\n0.000000e+00 0.000000e+00 "
4288  << std::sqrt(2/3.)/pow(det, 1./3.)
4289  << std::endl;
4290 
4291  if (me == 3)
4292  metric_file << 1./g1 << " "
4293  << 0.5/g2 << " "
4294  << 0.5/g3 << "\n0.000000e+00 "
4295  << 0.5*std::sqrt(3.)/g2 << " "
4296  << 0.5/(std::sqrt(3.)*g3) << "\n0.000000e+00 0.000000e+00 "
4297  << std::sqrt(2./3.)/g3
4298  << std::endl;
4299 
4300  det_o = det;
4301  g1_o = g1;
4302  g2_o = g2;
4303  g3_o = g3;
4304  Ncells++;
4305  } // end for ni
4306  } // end hex
4307  } // end dim==3
4308  }
4309 
4310  // Done with the metric file
4311  metric_file.close();
4312 
4313  std::ofstream grid_file(grid.c_str());
4314  if (!grid_file.good())
4315  libmesh_error_msg("Error opening file: " << grid);
4316 
4317  grid_file << _dim << "\n" << _n_nodes << "\n" << Ncells << "\n0" << std::endl;
4318 
4319  // Use scientific notation with 6 digits
4320  grid_file << std::scientific << std::setprecision(6);
4321 
4322  for (dof_id_type i=0; i<_n_nodes; i++)
4323  {
4324  // node coordinates
4325  for (unsigned j=0; j<_dim; j++)
4326  grid_file << R[i][j] << " ";
4327 
4328  grid_file << mask[i] << std::endl;
4329  }
4330 
4331  for (dof_id_type i=0; i<_n_cells; i++)
4332  if (mcells[i] >= 0)
4333  {
4334  // cell connectivity
4335  int nvert = 0;
4336  while (cells[i][nvert] >= 0)
4337  nvert++;
4338 
4339  if ((nvert == 3) || ((_dim == 3) && (nvert == 4)))
4340  {
4341  // tri & tetr
4342  for (int j=0; j<nvert; j++)
4343  grid_file << cells[i][j] << " ";
4344 
4345  for (unsigned j=nvert; j<3*_dim + _dim%2; j++)
4346  grid_file << "-1 ";
4347 
4348  grid_file << mcells[i] << std::endl;
4349  }
4350 
4351  if ((_dim == 2) && (nvert == 4))
4352  {
4353  // quad
4354  grid_file << cells[i][0] << " " << cells[i][1] << " " << cells[i][2] << " -1 -1 -1 0" << std::endl;
4355  grid_file << cells[i][1] << " " << cells[i][3] << " " << cells[i][0] << " -1 -1 -1 0" << std::endl;
4356  grid_file << cells[i][3] << " " << cells[i][2] << " " << cells[i][1] << " -1 -1 -1 0" << std::endl;
4357  grid_file << cells[i][2] << " " << cells[i][0] << " " << cells[i][3] << " -1 -1 -1 0" << std::endl;
4358  }
4359 
4360  if (nvert == 8)
4361  {
4362  // hex
4363  grid_file << cells[i][0] << " " << cells[i][1] << " " << cells[i][2] << " " << cells[i][4] << " -1 -1 -1 -1 -1 -1 0" << std::endl;
4364  grid_file << cells[i][1] << " " << cells[i][3] << " " << cells[i][0] << " " << cells[i][5] << " -1 -1 -1 -1 -1 -1 0" << std::endl;
4365  grid_file << cells[i][3] << " " << cells[i][2] << " " << cells[i][1] << " " << cells[i][7] << " -1 -1 -1 -1 -1 -1 0" << std::endl;
4366  grid_file << cells[i][2] << " " << cells[i][0] << " " << cells[i][3] << " " << cells[i][6] << " -1 -1 -1 -1 -1 -1 0" << std::endl;
4367  grid_file << cells[i][4] << " " << cells[i][6] << " " << cells[i][5] << " " << cells[i][0] << " -1 -1 -1 -1 -1 -1 0" << std::endl;
4368  grid_file << cells[i][5] << " " << cells[i][4] << " " << cells[i][7] << " " << cells[i][1] << " -1 -1 -1 -1 -1 -1 0" << std::endl;
4369  grid_file << cells[i][7] << " " << cells[i][5] << " " << cells[i][6] << " " << cells[i][3] << " -1 -1 -1 -1 -1 -1 0" << std::endl;
4370  grid_file << cells[i][6] << " " << cells[i][7] << " " << cells[i][4] << " " << cells[i][2] << " -1 -1 -1 -1 -1 -1 0" << std::endl;
4371  }
4372  }
4373 }
4374 
4375 } // namespace libMesh
4376 
4377 #endif // LIBMESH_ENABLE_VSMOOTHER
libMesh::Quality::name
std::string name(const ElemQuality q)
Definition: elem_quality.C:42
std::abs
double abs(double a)
Definition: libmesh_augment_std_namespace.h:59
libMesh::VariationalMeshSmoother::solver
int solver(int n, const std::vector< int > &ia, const std::vector< int > &ja, const std::vector< double > &a, std::vector< double > &x, const std::vector< double > &b, double eps, int maxite, int msglev)
Definition: mesh_smoother_vsmoother.C:3627
libMesh::VariationalMeshSmoother::_hanging_nodes
std::map< dof_id_type, std::vector< dof_id_type > > _hanging_nodes
Definition: mesh_smoother_vsmoother.h:169
libMesh::VariationalMeshSmoother::basisA
int basisA(Array2D< double > &Q, int nvert, const std::vector< double > &K, const Array2D< double > &H, int me)
Definition: mesh_smoother_vsmoother.C:621
libMesh::VariationalMeshSmoother::full_smooth
void full_smooth(Array2D< double > &R, const std::vector< int > &mask, const Array2D< int > &cells, const std::vector< int > &mcells, const std::vector< int > &edges, const std::vector< int > &hnodes, double w, const std::vector< int > &iter, int me, const Array3D< double > &H, int adp, int gr)
Definition: mesh_smoother_vsmoother.C:876
libMesh::VariationalMeshSmoother::adp_renew
void adp_renew(const Array2D< double > &R, const Array2D< int > &cells, std::vector< double > &afun, int adp)
Definition: mesh_smoother_vsmoother.C:829
libMesh::MeshBase::n_active_elem
virtual dof_id_type n_active_elem() const =0
libMesh::Node
A geometric point in (x,y,z) space associated with a DOF.
Definition: node.h:52
libMesh::VariationalMeshSmoother::jac2
double jac2(double x1, double y1, double x2, double y2)
Definition: mesh_smoother_vsmoother.C:610
libMesh::VariationalMeshSmoother::minJ
double minJ(Array2D< double > &R, const std::vector< int > &mask, const Array2D< int > &cells, const std::vector< int > &mcells, double epsilon, double w, int me, const Array3D< double > &H, double vol, const std::vector< int > &edges, const std::vector< int > &hnodes, int msglev, double &Vmin, double &emax, double &qmin, int adp, const std::vector< double > &afun)
Definition: mesh_smoother_vsmoother.C:1887
libMesh::MeshBase::const_element_iterator
Definition: mesh_base.h:1527
libMesh::VariationalMeshSmoother::vertex
double vertex(Array3D< double > &W, Array2D< double > &F, const Array2D< double > &R, const std::vector< int > &cell_in, double epsilon, double w, int nvert, const std::vector< double > &K, const Array2D< double > &H, int me, double vol, int f, double &Vmin, int adp, const std::vector< double > &g, double sigma)
Definition: mesh_smoother_vsmoother.C:3409
mesh
MeshBase & mesh
Definition: mesh_communication.C:1245
libMesh::MeshTools::find_nodal_neighbors
void find_nodal_neighbors(const MeshBase &mesh, const Node &n, const std::vector< std::vector< const Elem *>> &nodes_to_elem_map, std::vector< const Node *> &neighbors)
Definition: mesh_tools.C:740
libMesh::MeshTools::build_nodes_to_elem_map
void build_nodes_to_elem_map(const MeshBase &mesh, std::vector< std::vector< dof_id_type >> &nodes_to_elem_map)
Definition: mesh_tools.C:245
end
IterBase * end
Definition: variant_filter_iterator.h:343
libMesh::MeshBase::active_element_ptr_range
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
libMesh::MeshSmoother::_mesh
UnstructuredMesh & _mesh
Definition: mesh_smoother.h:61
std::max
long double max(long double a, double b)
Definition: libmesh_augment_std_namespace.h:27
libMesh::VariationalMeshSmoother::minq
double minq(const Array2D< double > &R, const Array2D< int > &cells, const std::vector< int > &mcells, int me, const Array3D< double > &H, double &vol, double &Vmin)
Definition: mesh_smoother_vsmoother.C:1484
libMesh::VariationalMeshSmoother::readgr
int readgr(Array2D< double > &R, std::vector< int > &mask, Array2D< int > &cells, std::vector< int > &mcells, std::vector< int > &edges, std::vector< int > &hnodes)
Definition: mesh_smoother_vsmoother.C:269
libMesh::VariationalMeshSmoother::minJ_BC
double minJ_BC(Array2D< double > &R, const std::vector< int > &mask, const Array2D< int > &cells, const std::vector< int > &mcells, double epsilon, double w, int me, const Array3D< double > &H, double vol, int msglev, double &Vmin, double &emax, double &qmin, int adp, const std::vector< double > &afun, int NCN)
Definition: mesh_smoother_vsmoother.C:2359
libMesh::MeshBase::elements_begin
virtual element_iterator elements_begin()=0
unstructured_mesh.h
libMesh::VariationalMeshSmoother::localP
double localP(Array3D< double > &W, Array2D< double > &F, Array2D< double > &R, const std::vector< int > &cell_in, const std::vector< int > &mask, double epsilon, double w, int nvert, const Array2D< double > &H, int me, double vol, int f, double &Vmin, double &qmin, int adp, const std::vector< double > &afun, std::vector< double > &Gloc)
Definition: mesh_smoother_vsmoother.C:2979
libMesh::UnstructuredMesh
Base class for Replicated and Distributed meshes.
Definition: unstructured_mesh.h:48
libMesh::MeshBase::elements_end
virtual element_iterator elements_end()=0
libMesh::MeshBase::node_ptr_range
virtual SimpleRange< node_iterator > node_ptr_range()=0
libMesh::VariationalMeshSmoother::pcg_par_check
int pcg_par_check(int n, const std::vector< int > &ia, const std::vector< int > &ja, const std::vector< double > &a, double eps, int maxite, int msglev)
Definition: mesh_smoother_vsmoother.C:3664
libMesh::VariationalMeshSmoother::jac3
double jac3(double x1, double y1, double z1, double x2, double y2, double z2, double x3, double y3, double z3)
Definition: mesh_smoother_vsmoother.C:595
libMesh::VariationalMeshSmoother::readmetr
int readmetr(std::string name, Array3D< double > &H)
Definition: mesh_smoother_vsmoother.C:490
libMesh::VariationalMeshSmoother::maxE
double maxE(Array2D< double > &R, const Array2D< int > &cells, const std::vector< int > &mcells, int me, const Array3D< double > &H, double v, double epsilon, double w, std::vector< double > &Gamma, double &qmin)
Definition: mesh_smoother_vsmoother.C:1077
libMesh::VariationalMeshSmoother::read_adp
int read_adp(std::vector< double > &afun)
Definition: mesh_smoother_vsmoother.C:573
libMesh::MeshTools::find_hanging_nodes_and_parents
void find_hanging_nodes_and_parents(const MeshBase &mesh, std::map< dof_id_type, std::vector< dof_id_type >> &hanging_nodes)
Definition: mesh_tools.C:1083
std::pow
double pow(double a, int b)
Definition: libmesh_augment_std_namespace.h:66
libMesh::VariationalMeshSmoother::VariationalMeshSmoother
VariationalMeshSmoother(UnstructuredMesh &mesh, double theta=0.5, unsigned miniter=2, unsigned maxiter=5, unsigned miniterBC=5)
Definition: mesh_smoother_vsmoother.C:45
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131
A
static PetscErrorCode Mat * A
Definition: petscdmlibmeshimpl.C:1020
libMesh::VariationalMeshSmoother::pcg_ic0
int pcg_ic0(int n, const std::vector< int > &ia, const std::vector< int > &ja, const std::vector< double > &a, const std::vector< double > &u, std::vector< double > &x, const std::vector< double > &b, std::vector< double > &r, std::vector< double > &p, std::vector< double > &z, double eps, int maxite, int msglev)
Definition: mesh_smoother_vsmoother.C:3806
mesh_tools.h
libmesh_config.h
libMesh::VariationalMeshSmoother::avertex
double avertex(const std::vector< double > &afun, std::vector< double > &G, const Array2D< double > &R, const std::vector< int > &cell_in, int nvert, int adp)
Definition: mesh_smoother_vsmoother.C:3317
libMesh::VariationalMeshSmoother::_area_of_interest
const UnstructuredMesh * _area_of_interest
Definition: mesh_smoother_vsmoother.h:217
libMesh::out
OStreamProxy out(std::cout)
Definition: libmesh_common.h:228
libMesh::VariationalMeshSmoother::writegr
int writegr(const Array2D< double > &R)
Definition: mesh_smoother_vsmoother.C:223
mesh_smoother_vsmoother.h
std::min
long double min(long double a, double b)
Definition: libmesh_augment_std_namespace.h:29
libMesh::Point
A geometric point in (x,y,z) space.
Definition: point.h:38
libMesh::MeshBase::n_nodes
virtual dof_id_type n_nodes() const =0
libMesh::MeshTools::find_boundary_nodes
void find_boundary_nodes(const MeshBase &mesh, std::vector< bool > &on_boundary)
Definition: mesh_tools.C:303
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:64
libMesh::pi
const Real pi
Definition: libmesh.h:233
libMesh::VariationalMeshSmoother::metr_data_gen
void metr_data_gen(std::string grid, std::string metr, int me)
Definition: mesh_smoother_vsmoother.C:3988
generated by doxygen