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