libMesh
Classes | Public Types | Public Member Functions | Protected Attributes | Private Member Functions | Private Attributes | List of all members
libMesh::VariationalMeshSmoother Class Reference

This is an implementation of Larisa Branets' smoothing algorithms. More...

#include <mesh_smoother_vsmoother.h>

Inheritance diagram for libMesh::VariationalMeshSmoother:
[legend]

Classes

struct  Array2D
 2D array type for interfacing with C APIs. More...
 
struct  Array3D
 3D array type for interfacing with C APIs. More...
 

Public Types

enum  MetricType { UNIFORM = 1, VOLUMETRIC = 2, DIRECTIONAL = 3 }
 
enum  AdaptType { CELL = -1, NONE = 0, NODE = 1 }
 

Public Member Functions

 VariationalMeshSmoother (UnstructuredMesh &mesh, Real theta=0.5, unsigned miniter=2, unsigned maxiter=5, unsigned miniterBC=5)
 Simple constructor to use for smoothing purposes. More...
 
 VariationalMeshSmoother (UnstructuredMesh &mesh, std::vector< float > *adapt_data, Real theta=0.5, unsigned miniter=2, unsigned maxiter=5, unsigned miniterBC=5, Real percent_to_move=1)
 Slightly more complicated constructor for mesh redistribution based on adapt_data. More...
 
 VariationalMeshSmoother (UnstructuredMesh &mesh, const UnstructuredMesh *area_of_interest, std::vector< float > *adapt_data, Real theta=0.5, unsigned miniter=2, unsigned maxiter=5, unsigned miniterBC=5, Real percent_to_move=1)
 Even more complicated constructor for mesh redistribution based on adapt_data with an area of interest. More...
 
virtual ~VariationalMeshSmoother ()
 Destructor. More...
 
virtual void smooth () override
 Redefinition of the smooth function from the base class. More...
 
Real smooth (unsigned int n_iterations)
 The actual smoothing function, gets called whenever the user specifies an actual number of smoothing iterations. More...
 
Real distance_moved () const
 
void set_generate_data (bool b)
 Allow user to control whether the metric is generated from the initial mesh. More...
 
void set_metric (MetricType t)
 Allow user to control the smoothing metric used. More...
 

Protected Attributes

UnstructuredMesh & _mesh
 

Private Member Functions

void adjust_adapt_data ()
 
float adapt_minimum () const
 
int writegr (const Array2D< Real > &R)
 
int readgr (Array2D< Real > &R, std::vector< int > &mask, Array2D< int > &cells, std::vector< int > &mcells, std::vector< int > &edges, std::vector< int > &hnodes)
 
int readmetr (std::string name, Array3D< Real > &H)
 
int read_adp (std::vector< Real > &afun)
 
Real jac3 (Real x1, Real y1, Real z1, Real x2, Real y2, Real z2, Real x3, Real y3, Real z3)
 
Real jac2 (Real x1, Real y1, Real x2, Real y2)
 
int basisA (Array2D< Real > &Q, int nvert, const std::vector< Real > &K, const Array2D< Real > &H, int me)
 
void adp_renew (const Array2D< Real > &R, const Array2D< int > &cells, std::vector< Real > &afun, int adp)
 
void full_smooth (Array2D< Real > &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, Real w, const std::vector< int > &iter, int me, const Array3D< Real > &H, int adp, int gr)
 
Real maxE (Array2D< Real > &R, const Array2D< int > &cells, const std::vector< int > &mcells, int me, const Array3D< Real > &H, Real v, Real epsilon, Real w, std::vector< Real > &Gamma, Real &qmin)
 
Real minq (const Array2D< Real > &R, const Array2D< int > &cells, const std::vector< int > &mcells, int me, const Array3D< Real > &H, Real &vol, Real &Vmin)
 
Real minJ (Array2D< Real > &R, const std::vector< int > &mask, const Array2D< int > &cells, const std::vector< int > &mcells, Real epsilon, Real w, int me, const Array3D< Real > &H, Real vol, const std::vector< int > &edges, const std::vector< int > &hnodes, int msglev, Real &Vmin, Real &emax, Real &qmin, int adp, const std::vector< Real > &afun)
 
Real minJ_BC (Array2D< Real > &R, const std::vector< int > &mask, const Array2D< int > &cells, const std::vector< int > &mcells, Real epsilon, Real w, int me, const Array3D< Real > &H, Real vol, int msglev, Real &Vmin, Real &emax, Real &qmin, int adp, const std::vector< Real > &afun, int NCN)
 
Real localP (Array3D< Real > &W, Array2D< Real > &F, Array2D< Real > &R, const std::vector< int > &cell_in, const std::vector< int > &mask, Real epsilon, Real w, int nvert, const Array2D< Real > &H, int me, Real vol, int f, Real &Vmin, Real &qmin, int adp, const std::vector< Real > &afun, std::vector< Real > &Gloc)
 
Real avertex (const std::vector< Real > &afun, std::vector< Real > &G, const Array2D< Real > &R, const std::vector< int > &cell_in, int nvert, int adp)
 
Real vertex (Array3D< Real > &W, Array2D< Real > &F, const Array2D< Real > &R, const std::vector< int > &cell_in, Real epsilon, Real w, int nvert, const std::vector< Real > &K, const Array2D< Real > &H, int me, Real vol, int f, Real &Vmin, int adp, const std::vector< Real > &g, Real sigma)
 
void metr_data_gen (std::string grid, std::string metr, int me)
 
int solver (int n, const std::vector< int > &ia, const std::vector< int > &ja, const std::vector< Real > &a, std::vector< Real > &x, const std::vector< Real > &b, Real eps, int maxite, int msglev)
 
int pcg_ic0 (int n, const std::vector< int > &ia, const std::vector< int > &ja, const std::vector< Real > &a, const std::vector< Real > &u, std::vector< Real > &x, const std::vector< Real > &b, std::vector< Real > &r, std::vector< Real > &p, std::vector< Real > &z, Real eps, int maxite, int msglev)
 
int pcg_par_check (int n, const std::vector< int > &ia, const std::vector< int > &ja, const std::vector< Real > &a, Real eps, int maxite, int msglev)
 
void gener (char grid[], int n)
 

Private Attributes

Real _distance
 Max distance of the last set of movement. More...
 
const Real _percent_to_move
 Dampening factor. More...
 
Real _dist_norm
 Records a relative "distance moved". More...
 
std::map< dof_id_type, std::vector< dof_id_type > > _hanging_nodes
 Map for hanging_nodes. More...
 
std::vector< float > * _adapt_data
 Vector for holding adaptive data. More...
 
const unsigned _dim
 Smoother control variables. More...
 
const unsigned _miniter
 
const unsigned _maxiter
 
const unsigned _miniterBC
 
MetricType _metric
 
AdaptType _adaptive_func
 
const Real _theta
 
bool _generate_data
 
dof_id_type _n_nodes
 The number of nodes in the Mesh at the time of smoothing. More...
 
dof_id_type _n_cells
 The number of active elements in the Mesh at the time of smoothing. More...
 
dof_id_type _n_hanging_edges
 The number of hanging node edges in the Mesh at the time of smoothing. More...
 
std::ofstream _logfile
 All output (including debugging) is sent to the _logfile. More...
 
const UnstructuredMesh * _area_of_interest
 Area of Interest Mesh. More...
 

Detailed Description

This is an implementation of Larisa Branets' smoothing algorithms.

The initial implementation was done by her, the adaptation to libmesh was completed by Derek Gaston. The code was heavily refactored into something more closely resembling C++ by John Peterson in 2014.

Here are the relevant publications: 1) L. Branets, G. Carey, "Extension of a mesh quality metric for elements with a curved boundary edge or surface", Journal of Computing and Information Science in Engineering, vol. 5(4), pp.302-308, 2005.

2) L. Branets, G. Carey, "A local cell quality metric and variational grid smoothing algorithm", Engineering with Computers, vol. 21, pp.19-28, 2005.

3) L. Branets, "A variational grid optimization algorithm based on a local cell quality metric", Ph.D. thesis, The University of Texas at Austin, 2005.

Author
Derek R. Gaston
Date
2006

Definition at line 63 of file mesh_smoother_vsmoother.h.

Member Enumeration Documentation

◆ AdaptType

enum libMesh::VariationalMeshSmoother::AdaptType
Enumerator
CELL 
NONE 
NODE 

Definition at line 107 of file mesh_smoother_vsmoother.h.

108  {
109  CELL = -1,
110  NONE = 0,
111  NODE = 1
112  };

◆ MetricType

enum libMesh::VariationalMeshSmoother::MetricType
Enumerator
UNIFORM 
VOLUMETRIC 
DIRECTIONAL 

Definition at line 100 of file mesh_smoother_vsmoother.h.

101  {
102  UNIFORM = 1,
103  VOLUMETRIC = 2,
104  DIRECTIONAL = 3
105  };

Constructor & Destructor Documentation

◆ VariationalMeshSmoother() [1/3]

libMesh::VariationalMeshSmoother::VariationalMeshSmoother ( UnstructuredMesh &  mesh,
Real  theta = 0.5,
unsigned  miniter = 2,
unsigned  maxiter = 5,
unsigned  miniterBC = 5 
)

Simple constructor to use for smoothing purposes.

Definition at line 45 of file mesh_smoother_vsmoother.C.

49  :
50  MeshSmoother(mesh),
51  _percent_to_move(1),
52  _dist_norm(0.),
53  _adapt_data(nullptr),
54  _dim(mesh.mesh_dimension()),
55  _miniter(miniter),
56  _maxiter(maxiter),
57  _miniterBC(miniterBC),
58  _metric(UNIFORM),
59  _adaptive_func(NONE),
60  _theta(theta),
61  _generate_data(false),
62  _n_nodes(0),
63  _n_cells(0),
64  _n_hanging_edges(0),
65  _area_of_interest(nullptr)
66 {}

◆ VariationalMeshSmoother() [2/3]

libMesh::VariationalMeshSmoother::VariationalMeshSmoother ( UnstructuredMesh &  mesh,
std::vector< float > *  adapt_data,
Real  theta = 0.5,
unsigned  miniter = 2,
unsigned  maxiter = 5,
unsigned  miniterBC = 5,
Real  percent_to_move = 1 
)

Slightly more complicated constructor for mesh redistribution based on adapt_data.

Definition at line 71 of file mesh_smoother_vsmoother.C.

77  :
78  MeshSmoother(mesh),
79  _percent_to_move(percent_to_move),
80  _dist_norm(0.),
81  _adapt_data(adapt_data),
82  _dim(mesh.mesh_dimension()),
83  _miniter(miniter),
84  _maxiter(maxiter),
85  _miniterBC(miniterBC),
86  _metric(UNIFORM),
87  _adaptive_func(CELL),
88  _theta(theta),
89  _generate_data(false),
90  _n_nodes(0),
91  _n_cells(0),
92  _n_hanging_edges(0),
93  _area_of_interest(nullptr)
94 {}

◆ VariationalMeshSmoother() [3/3]

libMesh::VariationalMeshSmoother::VariationalMeshSmoother ( UnstructuredMesh &  mesh,
const UnstructuredMesh *  area_of_interest,
std::vector< float > *  adapt_data,
Real  theta = 0.5,
unsigned  miniter = 2,
unsigned  maxiter = 5,
unsigned  miniterBC = 5,
Real  percent_to_move = 1 
)

Even more complicated constructor for mesh redistribution based on adapt_data with an area of interest.

Definition at line 98 of file mesh_smoother_vsmoother.C.

105  :
106  MeshSmoother(mesh),
107  _percent_to_move(percent_to_move),
108  _dist_norm(0.),
109  _adapt_data(adapt_data),
110  _dim(mesh.mesh_dimension()),
111  _miniter(miniter),
112  _maxiter(maxiter),
113  _miniterBC(miniterBC),
114  _metric(UNIFORM),
115  _adaptive_func(CELL),
116  _theta(theta),
117  _generate_data(false),
118  _n_nodes(0),
119  _n_cells(0),
120  _n_hanging_edges(0),
121  _area_of_interest(area_of_interest)
122 {}

◆ ~VariationalMeshSmoother()

virtual libMesh::VariationalMeshSmoother::~VariationalMeshSmoother ( )
inlinevirtual

Destructor.

Definition at line 117 of file mesh_smoother_vsmoother.h.

117 {}

Member Function Documentation

◆ adapt_minimum()

float libMesh::VariationalMeshSmoother::adapt_minimum ( ) const
private

Definition at line 512 of file mesh_smoother_vsmoother.C.

513 {
514  float min = std::numeric_limits<float>::max();
515 
516  for (const auto & val : *_adapt_data)
517  {
518  // Only positive (or zero) values in the error vector
519  libmesh_assert_greater_equal (val, 0.);
520  min = std::min (min, val);
521  }
522 
523  // ErrorVectors are for positive values
524  libmesh_assert_greater_equal (min, 0.);
525 
526  return min;
527 }

References _adapt_data.

Referenced by adjust_adapt_data().

◆ adjust_adapt_data()

void libMesh::VariationalMeshSmoother::adjust_adapt_data ( )
private

Definition at line 531 of file mesh_smoother_vsmoother.C.

532 {
533  // For convenience
534  const UnstructuredMesh & aoe_mesh = *_area_of_interest;
535  std::vector<float> & adapt_data = *_adapt_data;
536 
537  float min = adapt_minimum();
538 
539  MeshBase::const_element_iterator el = _mesh.elements_begin();
540  const MeshBase::const_element_iterator end_el = _mesh.elements_end();
541 
542  MeshBase::const_element_iterator aoe_el = aoe_mesh.elements_begin();
543  const MeshBase::const_element_iterator aoe_end_el = aoe_mesh.elements_end();
544 
545  // Counter to keep track of which spot in adapt_data we're looking at
546  for (int i=0; el!=end_el; el++)
547  {
548  // Only do this for active elements
549  if (adapt_data[i])
550  {
551  Point centroid = (*el)->centroid();
552  bool in_aoe = false;
553 
554  // See if the elements centroid lies in the aoe mesh
555  // for (aoe_el=aoe_mesh.elements_begin(); aoe_el != aoe_end_el; ++aoe_el)
556  // {
557  if ((*aoe_el)->contains_point(centroid))
558  in_aoe = true;
559  // }
560 
561  // If the element is not in the area of interest... then set
562  // it's adaptivity value to the minimum.
563  if (!in_aoe)
564  adapt_data[i] = min;
565  }
566  i++;
567  }
568 }

References _adapt_data, _area_of_interest, libMesh::MeshSmoother::_mesh, adapt_minimum(), libMesh::MeshBase::elements_begin(), and libMesh::MeshBase::elements_end().

Referenced by read_adp().

◆ adp_renew()

void libMesh::VariationalMeshSmoother::adp_renew ( const Array2D< Real > &  R,
const Array2D< int > &  cells,
std::vector< Real > &  afun,
int  adp 
)
private

Definition at line 829 of file mesh_smoother_vsmoother.C.

833 {
834  // evaluates given adaptive function on the provided mesh
835  if (adp < 0)
836  {
837  for (dof_id_type i=0; i<_n_cells; i++)
838  {
839  Real
840  x = 0.,
841  y = 0.,
842  z = 0.;
843  int nvert = 0;
844 
845  while (cells[i][nvert] >= 0)
846  {
847  x += R[cells[i][nvert]][0];
848  y += R[cells[i][nvert]][1];
849  if (_dim == 3)
850  z += R[cells[i][nvert]][2];
851  nvert++;
852  }
853 
854  // adaptive function, cell based
855  afun[i] = 5.*(x+y+z);
856  }
857  }
858 
859  else if (adp > 0)
860  {
861  for (dof_id_type i=0; i<_n_nodes; i++)
862  {
863  Real z = 0;
864  for (unsigned j=0; j<_dim; j++)
865  z += R[i][j];
866 
867  // adaptive function, node based
868  afun[i] = 5*std::sin(R[i][0]);
869  }
870  }
871 }

References _dim, _n_cells, _n_nodes, and libMesh::Real.

Referenced by full_smooth().

◆ avertex()

Real libMesh::VariationalMeshSmoother::avertex ( const std::vector< Real > &  afun,
std::vector< Real > &  G,
const Array2D< Real > &  R,
const std::vector< int > &  cell_in,
int  nvert,
int  adp 
)
private

Definition at line 3317 of file mesh_smoother_vsmoother.C.

3323 {
3324  std::vector<Real> a1(3), a2(3), a3(3), qu(3), K(8);
3325  Array2D<Real> Q(3, nvert);
3326 
3327  for (int i=0; i<8; i++)
3328  K[i] = 0.5; // cell center
3329 
3330  basisA(Q, nvert, K, Q, 1);
3331 
3332  for (unsigned i=0; i<_dim; i++)
3333  for (int j=0; j<nvert; j++)
3334  {
3335  a1[i] += Q[i][j]*R[cell_in[j]][0];
3336  a2[i] += Q[i][j]*R[cell_in[j]][1];
3337  if (_dim == 3)
3338  a3[i] += Q[i][j]*R[cell_in[j]][2];
3339 
3340  qu[i] += Q[i][j]*afun[cell_in[j]];
3341  }
3342 
3343  Real det = 0.;
3344 
3345  if (_dim == 2)
3346  det = jac2(a1[0], a1[1],
3347  a2[0], a2[1]);
3348  else
3349  det = jac3(a1[0], a1[1], a1[2],
3350  a2[0], a2[1], a2[2],
3351  a3[0], a3[1], a3[2]);
3352 
3353  // Return val
3354  Real g = 0.;
3355 
3356  if (det != 0)
3357  {
3358  if (_dim == 2)
3359  {
3360  Real df0 = jac2(qu[0], qu[1], a2[0], a2[1])/det;
3361  Real df1 = jac2(a1[0], a1[1], qu[0], qu[1])/det;
3362  g = (1. + df0*df0 + df1*df1);
3363 
3364  if (adp == 2)
3365  {
3366  // directional adaptation G=diag(g_i)
3367  G[0] = 1. + df0*df0;
3368  G[1] = 1. + df1*df1;
3369  }
3370  else
3371  {
3372  for (unsigned i=0; i<_dim; i++)
3373  G[i] = g; // simple adaptation G=gI
3374  }
3375  }
3376  else
3377  {
3378  Real df0 = (qu[0]*(a2[1]*a3[2]-a2[2]*a3[1]) + qu[1]*(a2[2]*a3[0]-a2[0]*a3[2]) + qu[2]*(a2[0]*a3[1]-a2[1]*a3[0]))/det;
3379  Real df1 = (qu[0]*(a3[1]*a1[2]-a3[2]*a1[1]) + qu[1]*(a3[2]*a1[0]-a3[0]*a1[2]) + qu[2]*(a3[0]*a1[1]-a3[1]*a1[0]))/det;
3380  Real df2 = (qu[0]*(a1[1]*a2[2]-a1[2]*a2[1]) + qu[1]*(a1[2]*a2[0]-a1[0]*a2[2]) + qu[2]*(a1[0]*a2[1]-a1[1]*a2[0]))/det;
3381  g = 1. + df0*df0 + df1*df1 + df2*df2;
3382  if (adp == 2)
3383  {
3384  // directional adaptation G=diag(g_i)
3385  G[0] = 1. + df0*df0;
3386  G[1] = 1. + df1*df1;
3387  G[2] = 1. + df2*df2;
3388  }
3389  else
3390  {
3391  for (unsigned i=0; i<_dim; i++)
3392  G[i] = g; // simple adaptation G=gI
3393  }
3394  }
3395  }
3396  else
3397  {
3398  g = 1.;
3399  for (unsigned i=0; i<_dim; i++)
3400  G[i] = g;
3401  }
3402 
3403  return g;
3404 }

References _dim, basisA(), jac2(), jac3(), and libMesh::Real.

Referenced by localP().

◆ basisA()

int libMesh::VariationalMeshSmoother::basisA ( Array2D< Real > &  Q,
int  nvert,
const std::vector< Real > &  K,
const Array2D< Real > &  H,
int  me 
)
private

Definition at line 621 of file mesh_smoother_vsmoother.C.

626 {
627  Array2D<Real> U(_dim, nvert);
628 
629  // Some useful constants
630  const Real
631  sqrt3 = std::sqrt(3.),
632  sqrt6 = std::sqrt(6.);
633 
634  if (_dim == 2)
635  {
636  if (nvert == 4)
637  {
638  // quad
639  U[0][0] = -(1-K[1]);
640  U[0][1] = (1-K[1]);
641  U[0][2] = -K[1];
642  U[0][3] = K[1];
643 
644  U[1][0] = -(1-K[0]);
645  U[1][1] = -K[0];
646  U[1][2] = (1-K[0]);
647  U[1][3] = K[0];
648  }
649  else if (nvert == 3)
650  {
651  // tri
652  // for right target triangle
653  // U[0][0] = -1.; U[1][0] = -1.;
654  // U[0][1] = 1.; U[1][1] = 0.;
655  // U[0][2] = 0.; U[1][2] = 1.;
656 
657  // for regular triangle
658  U[0][0] = -1.;
659  U[0][1] = 1.;
660  U[0][2] = 0;
661 
662  U[1][0] = -1./sqrt3;
663  U[1][1] = -1./sqrt3;
664  U[1][2] = 2./sqrt3;
665  }
666  else if (nvert == 6)
667  {
668  // curved triangle
669  U[0][0] = -3*(1-K[0]-K[1]*0.5)+(K[0]-0.5*K[1])+K[1];
670  U[0][1] = -(1-K[0]-K[1]*0.5)+(K[0]-0.5*K[1])*3-K[1];
671  U[0][2] = 0;
672  U[0][3] = K[1]*4;
673  U[0][4] = -4*K[1];
674  U[0][5] = 4*(1-K[0]-K[1]*0.5)-4*(K[0]-0.5*K[1]);
675 
676  U[1][0] = (1-K[2]-K[3]*0.5)*(-1.5)+(K[2]-0.5*K[3])*(0.5)+K[3]*(0.5);
677  U[1][1] = (1-K[2]-K[3]*0.5)*(0.5)+(K[2]-0.5*K[3])*(-1.5)+K[3]*(0.5);
678  U[1][2] = (1-K[2]-K[3]*0.5)*(-1)+(K[2]-0.5*K[3])*(-1)+K[3]*(3);
679  U[1][3] = (K[2]-0.5*K[3])*(4.)+K[3]*(-2.);
680  U[1][4] = (1-K[2]-K[3]*0.5)*(4.)+K[3]*(-2.);
681  U[1][5] = (1-K[2]-K[3]*0.5)*(-2.)+(K[2]-0.5*K[3])*(-2.);
682  }
683  else
684  libmesh_error_msg("Invalid nvert = " << nvert);
685  }
686 
687  if (_dim == 3)
688  {
689  if (nvert == 8)
690  {
691  // hex
692  U[0][0] = -(1-K[0])*(1-K[1]);
693  U[0][1] = (1-K[0])*(1-K[1]);
694  U[0][2] = -K[0]*(1-K[1]);
695  U[0][3] = K[0]*(1-K[1]);
696  U[0][4] = -(1-K[0])*K[1];
697  U[0][5] = (1-K[0])*K[1];
698  U[0][6] = -K[0]*K[1];
699  U[0][7] = K[0]*K[1];
700 
701  U[1][0] = -(1-K[2])*(1-K[3]);
702  U[1][1] = -K[2]*(1-K[3]);
703  U[1][2] = (1-K[2])*(1-K[3]);
704  U[1][3] = K[2]*(1-K[3]);
705  U[1][4] = -(1-K[2])*K[3];
706  U[1][5] = -K[2]*K[3];
707  U[1][6] = (1-K[2])*K[3];
708  U[1][7] = K[2]*K[3];
709 
710  U[2][0] = -(1-K[4])*(1-K[5]);
711  U[2][1] = -K[4]*(1-K[5]);
712  U[2][2] = -(1-K[4])*K[5];
713  U[2][3] = -K[4]*K[5];
714  U[2][4] = (1-K[4])*(1-K[5]);
715  U[2][5] = K[4]*(1-K[5]);
716  U[2][6] = (1-K[4])*K[5];
717  U[2][7] = K[4]*K[5];
718  }
719  else if (nvert == 4)
720  {
721  // linear tetr
722  // for regular reference tetrahedron
723  U[0][0] = -1;
724  U[0][1] = 1;
725  U[0][2] = 0;
726  U[0][3] = 0;
727 
728  U[1][0] = -1./sqrt3;
729  U[1][1] = -1./sqrt3;
730  U[1][2] = 2./sqrt3;
731  U[1][3] = 0;
732 
733  U[2][0] = -1./sqrt6;
734  U[2][1] = -1./sqrt6;
735  U[2][2] = -1./sqrt6;
736  U[2][3] = 3./sqrt6;
737 
738  // for right corner reference tetrahedron
739  // U[0][0] = -1; U[1][0] = -1; U[2][0] = -1;
740  // U[0][1] = 1; U[1][1] = 0; U[2][1] = 0;
741  // U[0][2] = 0; U[1][2] = 1; U[2][2] = 0;
742  // U[0][3] = 0; U[1][3] = 0; U[2][3] = 1;
743  }
744  else if (nvert == 6)
745  {
746  // prism
747  // for regular triangle in the prism base
748  U[0][0] = -1+K[0];
749  U[0][1] = 1-K[0];
750  U[0][2] = 0;
751  U[0][3] = -K[0];
752  U[0][4] = K[0];
753  U[0][5] = 0;
754 
755  U[1][0] = 0.5*(-1+K[1]);
756  U[1][1] = 0.5*(-1+K[1]);
757  U[1][2] = (1-K[1]);
758  U[1][3] = -0.5*K[1];
759  U[1][4] = -0.5*K[1];
760  U[1][5] = K[1];
761 
762  U[2][0] = -1. + K[2] + 0.5*K[3];
763  U[2][1] = -K[2] + 0.5*K[3];
764  U[2][2] = -K[3];
765  U[2][3] = 1 - K[2] - 0.5*K[3];
766  U[2][4] = K[2] - 0.5*K[3];
767  U[2][5] = K[3];
768  }
769  else if (nvert == 10)
770  {
771  // quad tetr
772  U[0][0] = -3.*(1 - K[0] - 0.5*K[1] - K[2]/3.) + (K[0] - 0.5*K[1] - K[2]/3.) + (K[1] - K[2]/3.) + K[2];
773  U[0][1] = -1.*(1 - K[0] - 0.5*K[1] - K[2]/3.) + 3.*(K[0] - 0.5*K[1] - K[2]/3.) - (K[1] - K[2]/3.) - K[2];
774  U[0][2] = 0.;
775  U[0][3] = 0.;
776  U[0][4] = 4.*(K[1] - K[2]/3.);
777  U[0][5] = -4.*(K[1] - K[2]/3.);
778  U[0][6] = 4.*(1. - K[0] - 0.5*K[1] - K[2]/3.) - 4.*(K[0] - 0.5*K[1] - K[2]/3.);
779  U[0][7] = 4.*K[2];
780  U[0][8] = 0.;
781  U[0][9] = -4.*K[2];
782 
783  U[1][0] = -1.5*(1. - K[3] - K[4]*0.5 - K[5]/3.) + 0.5*(K[3] - 0.5*K[4] - K[5]/3.) + 0.5*(K[4] - K[5]/3.) + 0.5*K[5];
784  U[1][1] = 0.5*(1. - K[3] - K[4]*0.5 - K[5]/3.) - 1.5*(K[3] - 0.5*K[4] - K[5]/3.) + 0.5*(K[4] - K[5]/3.) + 0.5*K[5];
785  U[1][2] = -1.*(1. - K[3] - K[4]*0.5 - K[5]/3.) - (K[3] - 0.5*K[4] - K[5]/3.) + 3.*(K[4] - K[5]/3.) - K[5];
786  U[1][3] = 0.;
787  U[1][4] = 4.*( K[3] - 0.5*K[4] - K[5]/3.) - 2.*(K[4] - K[5]/3.);
788  U[1][5] = 4.*(1. - K[3] - 0.5*K[4] - K[5]/3.) - 2.*(K[4] - K[5]/3.);
789  U[1][6] = -2.*(1. - K[3] - 0.5*K[4] - K[5]/3.) - 2.*(K[3] - 0.5*K[4] - K[5]/3.);
790  U[1][7] = -2.*K[5];
791  U[1][8] = 4.*K[5];
792  U[1][9] = -2.*K[5];
793 
794  U[2][0] = -(1. - K[6] - 0.5*K[7] - K[8]/3.) + (K[6] - 0.5*K[7] - K[8]/3.)/3. + (K[7] - K[8]/3.)/3. + K[8]/3.;
795  U[2][1] = (1. - K[6] - 0.5*K[7] - K[8]/3.)/3. - (K[6] - 0.5*K[7] - K[8]/3.) + (K[7] - K[8]/3.)/3. + K[8]/3.;
796  U[2][2] = (1. - K[6] - 0.5*K[7] - K[8]/3.)/3. + (K[6] - 0.5*K[7] - K[8]/3.)/3. - (K[7] - K[8]/3.) + K[8]/3.;
797  U[2][3] = -(1. - K[6] - 0.5*K[7] - K[8]/3.) - (K[6] - 0.5*K[7] - K[8]/3.) - (K[7] - K[8]/3.) + 3.*K[8];
798  U[2][4] = -4.*(K[6] - K[7]*0.5 - K[8]/3.)/3. - 4.*(K[7] - K[8]/3.)/3.;
799  U[2][5] = -4.*(1. - K[6] - K[7]*0.5 - K[8]/3.)/3. - 4.*( K[7] - K[8]/3.)/3.;
800  U[2][6] = -4.*(1. - K[6] - K[7]*0.5 - K[8]/3.)/3. - 4.*(K[6] - 0.5*K[7] - K[8]/3.)/3.;
801  U[2][7] = 4.*(K[6] - K[7]*0.5 - K[8]/3.) - 4.*K[8]/3.;
802  U[2][8] = 4.*(K[7] - K[8]/3.) - 4.*K[8]/3.;
803  U[2][9] = 4.*(1. - K[6] - K[7]*0.5 - K[8]/3.) - 4.*K[8]/3.;
804  }
805  else
806  libmesh_error_msg("Invalid nvert = " << nvert);
807  }
808 
809  if (me == 1)
810  Q = U;
811 
812  else
813  {
814  for (unsigned i=0; i<_dim; i++)
815  for (int j=0; j<nvert; j++)
816  {
817  Q[i][j] = 0;
818  for (unsigned k=0; k<_dim; k++)
819  Q[i][j] += H[k][i]*U[k][j];
820  }
821  }
822 
823  return 0;
824 }

References _dim, libMesh::Real, and std::sqrt().

Referenced by avertex(), maxE(), metr_data_gen(), minq(), and vertex().

◆ distance_moved()

Real libMesh::VariationalMeshSmoother::distance_moved ( ) const
inline
Returns
Max distance a node moved during the last smooth.

Definition at line 137 of file mesh_smoother_vsmoother.h.

137 { return _distance; }

References _distance.

◆ full_smooth()

void libMesh::VariationalMeshSmoother::full_smooth ( Array2D< Real > &  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,
Real  w,
const std::vector< int > &  iter,
int  me,
const Array3D< Real > &  H,
int  adp,
int  gr 
)
private

Definition at line 876 of file mesh_smoother_vsmoother.C.

888 {
889  // Control the amount of print statements in this function
890  int msglev = 1;
891 
892  dof_id_type afun_size = 0;
893 
894  // Adaptive function is on cells
895  if (adp < 0)
896  afun_size = _n_cells;
897 
898  // Adaptive function is on nodes
899  else if (adp > 0)
900  afun_size = _n_nodes;
901 
902  std::vector<Real> afun(afun_size);
903  std::vector<int> maskf(_n_nodes);
904  std::vector<Real> Gamma(_n_cells);
905 
906  if (msglev >= 1)
907  _logfile << "N=" << _n_nodes
908  << " ncells=" << _n_cells
909  << " nedges=" << _n_hanging_edges
910  << std::endl;
911 
912 
913  // Boundary node counting
914  int NBN=0;
915  for (dof_id_type i=0; i<_n_nodes; i++)
916  if (mask[i] == 2 || mask[i] == 1)
917  NBN++;
918 
919  if (NBN > 0)
920  {
921  if (msglev >= 1)
922  _logfile << "# of Boundary Nodes=" << NBN << std::endl;
923 
924  NBN=0;
925  for (dof_id_type i=0; i<_n_nodes; i++)
926  if (mask[i] == 2)
927  NBN++;
928 
929  if (msglev >= 1)
930  _logfile << "# of moving Boundary Nodes=" << NBN << std::endl;
931  }
932 
933  for (dof_id_type i=0; i<_n_nodes; i++)
934  {
935  if ((mask[i]==1) || (mask[i]==2))
936  maskf[i] = 1;
937  else
938  maskf[i] = 0;
939  }
940 
941  // determination of min jacobian
942  Real vol, Vmin;
943  Real qmin = minq(R, cells, mcells, me, H, vol, Vmin);
944 
945  if (me > 1)
946  vol = 1.;
947 
948  if (msglev >= 1)
949  _logfile << "vol=" << vol
950  << " qmin=" << qmin
951  << " min volume = " << Vmin
952  << std::endl;
953 
954  // compute max distortion measure over all cells
955  Real epsilon = 1.e-9;
956  Real eps = qmin < 0 ? std::sqrt(epsilon*epsilon+0.004*qmin*qmin*vol*vol) : epsilon;
957  Real emax = maxE(R, cells, mcells, me, H, vol, eps, w, Gamma, qmin);
958 
959  if (msglev >= 1)
960  _logfile << " emax=" << emax << std::endl;
961 
962  // unfolding/smoothing
963 
964  // iteration tolerance
965  Real Enm1 = 1.;
966 
967  // read adaptive function from file
968  if (adp*gr != 0)
969  read_adp(afun);
970 
971  {
972  int
973  counter = 0,
974  ii = 0;
975 
976  while (((qmin <= 0) || (counter < iter[0]) || (std::abs(emax-Enm1) > 1e-3)) &&
977  (ii < iter[1]) &&
978  (counter < iter[1]))
979  {
980  libmesh_assert_less (counter, iter[1]);
981 
982  Enm1 = emax;
983 
984  if ((ii >= 0) && (ii % 2 == 0))
985  {
986  if (qmin < 0)
987  eps = std::sqrt(epsilon*epsilon + 0.004*qmin*qmin*vol*vol);
988  else
989  eps = epsilon;
990  }
991 
992  int ladp = adp;
993 
994  if ((qmin <= 0) || (counter < ii))
995  ladp = 0;
996 
997  // update adaptation function before each iteration
998  if ((ladp != 0) && (gr == 0))
999  adp_renew(R, cells, afun, adp);
1000 
1001  Real Jk = minJ(R, maskf, cells, mcells, eps, w, me, H, vol, edges, hnodes,
1002  msglev, Vmin, emax, qmin, ladp, afun);
1003 
1004  if (qmin > 0)
1005  counter++;
1006  else
1007  ii++;
1008 
1009  if (msglev >= 1)
1010  _logfile << "niter=" << counter
1011  << ", qmin*G/vol=" << qmin
1012  << ", Vmin=" << Vmin
1013  << ", emax=" << emax
1014  << ", Jk=" << Jk
1015  << ", Enm1=" << Enm1
1016  << std::endl;
1017  }
1018  }
1019 
1020  // BN correction - 2D only!
1021  epsilon = 1.e-9;
1022  if (NBN > 0)
1023  for (int counter=0; counter<iter[2]; counter++)
1024  {
1025  // update adaptation function before each iteration
1026  if ((adp != 0) && (gr == 0))
1027  adp_renew(R, cells, afun, adp);
1028 
1029  Real Jk = minJ_BC(R, mask, cells, mcells, eps, w, me, H, vol, msglev, Vmin, emax, qmin, adp, afun, NBN);
1030 
1031  if (msglev >= 1)
1032  _logfile << "NBC niter=" << counter
1033  << ", qmin*G/vol=" << qmin
1034  << ", Vmin=" << Vmin
1035  << ", emax=" << emax
1036  << std::endl;
1037 
1038  // Outrageous Enm1 to make sure we hit this at least once
1039  Enm1 = 99999;
1040 
1041  // Now that we've moved the boundary nodes (or not) we need to resmooth
1042  for (int j=0; j<iter[1]; j++)
1043  {
1044  if (std::abs(emax-Enm1) < 1e-2)
1045  break;
1046 
1047  // Save off the error from the previous smoothing step
1048  Enm1 = emax;
1049 
1050  // update adaptation function before each iteration
1051  if ((adp != 0) && (gr == 0))
1052  adp_renew(R, cells, afun, adp);
1053 
1054  Jk = minJ(R, maskf, cells, mcells, eps, w, me, H, vol, edges, hnodes, msglev, Vmin, emax, qmin, adp, afun);
1055 
1056  if (msglev >= 1)
1057  _logfile << " Re-smooth: niter=" << j
1058  << ", qmin*G/vol=" << qmin
1059  << ", Vmin=" << Vmin
1060  << ", emax=" << emax
1061  << ", Jk=" << Jk
1062  << std::endl;
1063  }
1064 
1065  if (msglev >= 1)
1066  _logfile << "NBC smoothed niter=" << counter
1067  << ", qmin*G/vol=" << qmin
1068  << ", Vmin=" << Vmin
1069  << ", emax=" << emax
1070  << std::endl;
1071  }
1072 }

References _logfile, _n_cells, _n_hanging_edges, _n_nodes, std::abs(), adp_renew(), maxE(), minJ(), minJ_BC(), minq(), read_adp(), libMesh::Real, and std::sqrt().

Referenced by smooth().

◆ gener()

void libMesh::VariationalMeshSmoother::gener ( char  grid[],
int  n 
)
private

Definition at line 3918 of file mesh_smoother_vsmoother.C.

3920 {
3921  const int n1 = 3;
3922 
3923  int
3924  N = 1,
3925  ncells = 1;
3926 
3927  for (int i=0; i<n; i++)
3928  {
3929  N *= n1;
3930  ncells *= (n1-1);
3931  }
3932 
3933  const Real x = 1./static_cast<Real>(n1-1);
3934 
3935  std::ofstream outfile(grid);
3936 
3937  outfile << n << "\n" << N << "\n" << ncells << "\n0\n" << std::endl;
3938 
3939  for (int i=0; i<N; i++)
3940  {
3941  // node coordinates
3942  int k = i;
3943  int mask = 0;
3944  for (int j=0; j<n; j++)
3945  {
3946  int ii = k % n1;
3947  if ((ii == 0) || (ii == n1-1))
3948  mask = 1;
3949 
3950  outfile << static_cast<Real>(ii)*x << " ";
3951  k /= n1;
3952  }
3953  outfile << mask << std::endl;
3954  }
3955 
3956  for (int i=0; i<ncells; i++)
3957  {
3958  // cell connectivity
3959  int nc = i;
3960  int ii = nc%(n1-1);
3961  nc /= (n1-1);
3962  int jj = nc%(n1-1);
3963  int kk = nc/(n1-1);
3964 
3965  if (n == 2)
3966  outfile << ii+n1*jj << " "
3967  << ii+1+jj*n1 << " "
3968  << ii+(jj+1)*n1 << " "
3969  << ii+1+(jj+1)*n1 << " ";
3970 
3971  if (n == 3)
3972  outfile << ii+n1*jj+n1*n1*kk << " "
3973  << ii+1+jj*n1+n1*n1*kk << " "
3974  << ii+(jj+1)*n1+n1*n1*kk << " "
3975  << ii+1+(jj+1)*n1+n1*n1*kk << " "
3976  << ii+n1*jj+n1*n1*(kk+1) << " "
3977  << ii+1+jj*n1+n1*n1*(kk+1) << " "
3978  << ii+(jj+1)*n1+n1*n1*(kk+1) << " "
3979  << ii+1+(jj+1)*n1+n1*n1*(kk+1) << " ";
3980 
3981  outfile << "-1 -1 0 \n";
3982  }
3983 }

References libMesh::Real.

◆ jac2()

Real libMesh::VariationalMeshSmoother::jac2 ( Real  x1,
Real  y1,
Real  x2,
Real  y2 
)
private

Definition at line 610 of file mesh_smoother_vsmoother.C.

614 {
615  return x1*y2 - x2*y1;
616 }

Referenced by avertex(), maxE(), metr_data_gen(), minq(), and vertex().

◆ jac3()

Real libMesh::VariationalMeshSmoother::jac3 ( Real  x1,
Real  y1,
Real  z1,
Real  x2,
Real  y2,
Real  z2,
Real  x3,
Real  y3,
Real  z3 
)
private

Definition at line 595 of file mesh_smoother_vsmoother.C.

604 {
605  return x1*(y2*z3 - y3*z2) + y1*(z2*x3 - z3*x2) + z1*(x2*y3 - x3*y2);
606 }

Referenced by avertex(), maxE(), metr_data_gen(), minq(), and vertex().

◆ localP()

Real libMesh::VariationalMeshSmoother::localP ( Array3D< Real > &  W,
Array2D< Real > &  F,
Array2D< Real > &  R,
const std::vector< int > &  cell_in,
const std::vector< int > &  mask,
Real  epsilon,
Real  w,
int  nvert,
const Array2D< Real > &  H,
int  me,
Real  vol,
int  f,
Real &  Vmin,
Real &  qmin,
int  adp,
const std::vector< Real > &  afun,
std::vector< Real > &  Gloc 
)
private

Definition at line 2979 of file mesh_smoother_vsmoother.C.

2996 {
2997  // K - determines approximation rule for local integral over the cell
2998  std::vector<Real> K(9);
2999 
3000  // f - flag, f=0 for determination of Hessian and gradient of the functional,
3001  // f=1 for determination of the functional value only;
3002  // adaptivity is determined on the first step for adp>0 (nodal based)
3003  if (f == 0)
3004  {
3005  if (adp > 0)
3006  avertex(afun, Gloc, R, cell_in, nvert, adp);
3007  if (adp == 0)
3008  {
3009  for (unsigned i=0; i<_dim; i++)
3010  Gloc[i] = 1.;
3011  }
3012  }
3013 
3014  Real
3015  sigma = 0.,
3016  fun = 0,
3017  gqmin = 1e32,
3018  g = 0, // Vmin
3019  lqmin = 0.;
3020 
3021  // cell integration depending on cell type
3022  if (_dim == 2)
3023  {
3024  // 2D
3025  if (nvert == 3)
3026  {
3027  // tri
3028  sigma = 1.;
3029  fun += vertex(W, F, R, cell_in, epsilon, w, nvert, K, H, me, vol, f, lqmin, adp, Gloc, sigma);
3030  g += sigma*lqmin;
3031  if (gqmin > lqmin)
3032  gqmin = lqmin;
3033  }
3034  if (nvert == 4)
3035  {
3036  //quad
3037  for (unsigned i=0; i<2; i++)
3038  {
3039  K[0] = i;
3040  for (unsigned j=0; j<2; j++)
3041  {
3042  K[1] = j;
3043  sigma = 0.25;
3044  fun += vertex(W, F, R, cell_in, epsilon, w, nvert, K, H, me,
3045  vol, f, lqmin, adp, Gloc, sigma);
3046  g += sigma*lqmin;
3047  if (gqmin > lqmin)
3048  gqmin = lqmin;
3049  }
3050  }
3051  }
3052  else
3053  {
3054  // quad tri
3055  for (unsigned i=0; i<3; i++)
3056  {
3057  K[0] = i*0.5;
3058  unsigned k = i/2;
3059  K[1] = static_cast<Real>(k);
3060 
3061  for (unsigned j=0; j<3; j++)
3062  {
3063  K[2] = j*0.5;
3064  k = j/2;
3065  K[3] = static_cast<Real>(k);
3066  if (i == j)
3067  sigma = 1./12.;
3068  else
3069  sigma = 1./24.;
3070 
3071  fun += vertex(W, F, R, cell_in, epsilon, w, nvert, K, H, me,
3072  vol, f, lqmin, adp, Gloc, sigma);
3073  g += sigma*lqmin;
3074  if (gqmin > lqmin)
3075  gqmin = lqmin;
3076  }
3077  }
3078  }
3079  }
3080  if (_dim == 3)
3081  {
3082  // 3D
3083  if (nvert == 4)
3084  {
3085  // tetr
3086  sigma = 1.;
3087  fun += vertex(W, F, R, cell_in, epsilon, w, nvert, K, H, me,
3088  vol, f, lqmin, adp, Gloc, sigma);
3089  g += sigma*lqmin;
3090  if (gqmin > lqmin)
3091  gqmin = lqmin;
3092  }
3093  if (nvert == 6)
3094  {
3095  //prism
3096  for (unsigned i=0; i<2; i++)
3097  {
3098  K[0] = i;
3099  for (unsigned j=0; j<2; j++)
3100  {
3101  K[1] = j;
3102  for (unsigned k=0; k<3; k++)
3103  {
3104  K[2] = 0.5*static_cast<Real>(k);
3105  K[3] = static_cast<Real>(k % 2);
3106  sigma = 1./12.;
3107  fun += vertex(W, F, R, cell_in, epsilon, w, nvert, K, H, me,
3108  vol, f, lqmin, adp, Gloc, sigma);
3109  g += sigma*lqmin;
3110  if (gqmin > lqmin)
3111  gqmin = lqmin;
3112  }
3113  }
3114  }
3115  }
3116  if (nvert == 8)
3117  {
3118  // hex
3119  for (unsigned i=0; i<2; i++)
3120  {
3121  K[0] = i;
3122  for (unsigned j=0; j<2; j++)
3123  {
3124  K[1] = j;
3125  for (unsigned k=0; k<2; k++)
3126  {
3127  K[2] = k;
3128  for (unsigned l=0; l<2; l++)
3129  {
3130  K[3] = l;
3131  for (unsigned m=0; m<2; m++)
3132  {
3133  K[4] = m;
3134  for (unsigned nn=0; nn<2; nn++)
3135  {
3136  K[5] = nn;
3137 
3138  if ((i==nn) && (j==l) && (k==m))
3139  sigma = 1./27.;
3140 
3141  if (((i==nn) && (j==l) && (k!=m)) ||
3142  ((i==nn) && (j!=l) && (k==m)) ||
3143  ((i!=nn) && (j==l) && (k==m)))
3144  sigma = 1./54.;
3145 
3146  if (((i==nn) && (j!=l) && (k!=m)) ||
3147  ((i!=nn) && (j!=l) && (k==m)) ||
3148  ((i!=nn) && (j==l) && (k!=m)))
3149  sigma = 1./108.;
3150 
3151  if ((i!=nn) && (j!=l) && (k!=m))
3152  sigma=1./216.;
3153 
3154  fun += vertex(W, F, R, cell_in, epsilon, w, nvert, K, H, me,
3155  vol, f, lqmin, adp, Gloc, sigma);
3156  g += sigma*lqmin;
3157 
3158  if (gqmin > lqmin)
3159  gqmin = lqmin;
3160  }
3161  }
3162  }
3163  }
3164  }
3165  }
3166  }
3167  else
3168  {
3169  // quad tetr
3170  for (unsigned i=0; i<4; i++)
3171  {
3172  for (unsigned j=0; j<4; j++)
3173  {
3174  for (unsigned k=0; k<4; k++)
3175  {
3176  switch (i)
3177  {
3178  case 0:
3179  K[0] = 0;
3180  K[1] = 0;
3181  K[2] = 0;
3182  break;
3183 
3184  case 1:
3185  K[0] = 1;
3186  K[1] = 0;
3187  K[2] = 0;
3188  break;
3189 
3190  case 2:
3191  K[0] = 0.5;
3192  K[1] = 1;
3193  K[2] = 0;
3194  break;
3195 
3196  case 3:
3197  K[0] = 0.5;
3198  K[1] = 1./3.;
3199  K[2] = 1;
3200  break;
3201 
3202  default:
3203  break;
3204  }
3205 
3206  switch (j)
3207  {
3208  case 0:
3209  K[3] = 0;
3210  K[4] = 0;
3211  K[5] = 0;
3212  break;
3213 
3214  case 1:
3215  K[3] = 1;
3216  K[4] = 0;
3217  K[5] = 0;
3218  break;
3219 
3220  case 2:
3221  K[3] = 0.5;
3222  K[4] = 1;
3223  K[5] = 0;
3224  break;
3225 
3226  case 3:
3227  K[3] = 0.5;
3228  K[4] = 1./3.;
3229  K[5] = 1;
3230  break;
3231 
3232  default:
3233  break;
3234  }
3235 
3236  switch (k)
3237  {
3238  case 0:
3239  K[6] = 0;
3240  K[7] = 0;
3241  K[8] = 0;
3242  break;
3243 
3244  case 1:
3245  K[6] = 1;
3246  K[7] = 0;
3247  K[8] = 0;
3248  break;
3249 
3250  case 2:
3251  K[6] = 0.5;
3252  K[7] = 1;
3253  K[8] = 0;
3254  break;
3255 
3256  case 3:
3257  K[6] = 0.5;
3258  K[7] = 1./3.;
3259  K[8] = 1;
3260  break;
3261 
3262  default:
3263  break;
3264  }
3265 
3266  if ((i==j) && (j==k))
3267  sigma = 1./120.;
3268 
3269  else if ((i==j) || (j==k) || (i==k))
3270  sigma = 1./360.;
3271 
3272  else
3273  sigma = 1./720.;
3274 
3275  fun += vertex(W, F, R, cell_in, epsilon, w, nvert, K, H, me,
3276  vol, f, lqmin, adp, Gloc, sigma);
3277 
3278  g += sigma*lqmin;
3279  if (gqmin > lqmin)
3280  gqmin = lqmin;
3281  }
3282  }
3283  }
3284  }
3285  }
3286 
3287  // fixed nodes correction
3288  for (int ii=0; ii<nvert; ii++)
3289  {
3290  if (mask[cell_in[ii]] == 1)
3291  {
3292  for (unsigned kk=0; kk<_dim; kk++)
3293  {
3294  for (int jj=0; jj<nvert; jj++)
3295  {
3296  W[kk][ii][jj] = 0;
3297  W[kk][jj][ii] = 0;
3298  }
3299 
3300  W[kk][ii][ii] = 1;
3301  F[kk][ii] = 0;
3302  }
3303  }
3304  }
3305  // end of fixed nodes correction
3306 
3307  // Set up return value references
3308  Vmin = g;
3309  qmin = gqmin/vol;
3310 
3311  return fun;
3312 }

References _dim, avertex(), libMesh::Real, and vertex().

Referenced by minJ(), and minJ_BC().

◆ maxE()

Real libMesh::VariationalMeshSmoother::maxE ( Array2D< Real > &  R,
const Array2D< int > &  cells,
const std::vector< int > &  mcells,
int  me,
const Array3D< Real > &  H,
Real  v,
Real  epsilon,
Real  w,
std::vector< Real > &  Gamma,
Real &  qmin 
)
private

Definition at line 1077 of file mesh_smoother_vsmoother.C.

1087 {
1088  Array2D<Real> Q(3, 3*_dim + _dim%2);
1089  std::vector<Real> K(9);
1090 
1091  Real
1092  gemax = -1.e32,
1093  vmin = 1.e32;
1094 
1095  for (dof_id_type ii=0; ii<_n_cells; ii++)
1096  if (mcells[ii] >= 0)
1097  {
1098  // Final value of E will be saved in Gamma at the end of this loop
1099  Real E = 0.;
1100 
1101  if (_dim == 2)
1102  {
1103  if (cells[ii][3] == -1)
1104  {
1105  // tri
1106  basisA(Q, 3, K, H[ii], me);
1107 
1108  std::vector<Real> a1(3), a2(3);
1109  for (int k=0; k<2; k++)
1110  for (int l=0; l<3; l++)
1111  {
1112  a1[k] += Q[k][l]*R[cells[ii][l]][0];
1113  a2[k] += Q[k][l]*R[cells[ii][l]][1];
1114  }
1115 
1116  Real det = jac2(a1[0], a1[1], a2[0], a2[1]);
1117  Real tr = 0.5*(a1[0]*a1[0] + a2[0]*a2[0] + a1[1]*a1[1] + a2[1]*a2[1]);
1118  Real chi = 0.5*(det+std::sqrt(det*det+epsilon*epsilon));
1119  E = (1-w)*tr/chi + 0.5*w*(v + det*det/v)/chi;
1120 
1121  if (E > gemax)
1122  gemax = E;
1123  if (vmin > det)
1124  vmin = det;
1125 
1126  }
1127  else if (cells[ii][4] == -1)
1128  {
1129  // quad
1130  for (int i=0; i<2; i++)
1131  {
1132  K[0] = i;
1133  for (int j=0; j<2; j++)
1134  {
1135  K[1] = j;
1136  basisA(Q, 4, K, H[ii], me);
1137 
1138  std::vector<Real> a1(3), a2(3);
1139  for (int k=0; k<2; k++)
1140  for (int l=0; l<4; l++)
1141  {
1142  a1[k] += Q[k][l]*R[cells[ii][l]][0];
1143  a2[k] += Q[k][l]*R[cells[ii][l]][1];
1144  }
1145 
1146  Real det = jac2(a1[0],a1[1],a2[0],a2[1]);
1147  Real tr = 0.5*(a1[0]*a1[0] + a2[0]*a2[0] + a1[1]*a1[1] + a2[1]*a2[1]);
1148  Real chi = 0.5*(det+std::sqrt(det*det+epsilon*epsilon));
1149  E += 0.25*((1-w)*tr/chi + 0.5*w*(v + det*det/v)/chi);
1150 
1151  if (vmin > det)
1152  vmin = det;
1153  }
1154  }
1155 
1156  if (E > gemax)
1157  gemax = E;
1158  }
1159  else
1160  {
1161  // quad tri
1162  for (int i=0; i<3; i++)
1163  {
1164  K[0] = i*0.5;
1165  int k = i/2;
1166  K[1] = static_cast<Real>(k);
1167 
1168  for (int j=0; j<3; j++)
1169  {
1170  K[2] = j*0.5;
1171  k = j/2;
1172  K[3] = static_cast<Real>(k);
1173 
1174  basisA(Q, 6, K, H[ii], me);
1175 
1176  std::vector<Real> a1(3), a2(3);
1177  for (int k2=0; k2<2; k2++)
1178  for (int l=0; l<6; l++)
1179  {
1180  a1[k2] += Q[k2][l]*R[cells[ii][l]][0];
1181  a2[k2] += Q[k2][l]*R[cells[ii][l]][1];
1182  }
1183 
1184  Real det = jac2(a1[0],a1[1],a2[0],a2[1]);
1185  Real sigma = 1./24.;
1186 
1187  if (i==j)
1188  sigma = 1./12.;
1189 
1190  Real tr = 0.5*(a1[0]*a1[0] + a2[0]*a2[0] + a1[1]*a1[1] + a2[1]*a2[1]);
1191  Real chi = 0.5*(det + std::sqrt(det*det + epsilon*epsilon));
1192  E += sigma*((1-w)*tr/chi + 0.5*w*(v + det*det/v)/chi);
1193  if (vmin > det)
1194  vmin = det;
1195  }
1196  }
1197 
1198  if (E > gemax)
1199  gemax = E;
1200  }
1201  }
1202 
1203  if (_dim == 3)
1204  {
1205  if (cells[ii][4] == -1)
1206  {
1207  // tetr
1208  basisA(Q, 4, K, H[ii], me);
1209 
1210  std::vector<Real> a1(3), a2(3), a3(3);
1211  for (int k=0; k<3; k++)
1212  for (int l=0; l<4; l++)
1213  {
1214  a1[k] += Q[k][l]*R[cells[ii][l]][0];
1215  a2[k] += Q[k][l]*R[cells[ii][l]][1];
1216  a3[k] += Q[k][l]*R[cells[ii][l]][2];
1217  }
1218 
1219  Real det = jac3(a1[0], a1[1], a1[2],
1220  a2[0], a2[1], a2[2],
1221  a3[0], a3[1], a3[2]);
1222  Real tr = 0.;
1223  for (int k=0; k<3; k++)
1224  tr += (a1[k]*a1[k] + a2[k]*a2[k] + a3[k]*a3[k])/3.;
1225 
1226  Real chi = 0.5*(det+std::sqrt(det*det+epsilon*epsilon));
1227  E = (1-w)*pow(tr,1.5)/chi + 0.5*w*(v + det*det/v)/chi;
1228 
1229  if (E > gemax)
1230  gemax = E;
1231 
1232  if (vmin > det)
1233  vmin = det;
1234  }
1235  else if (cells[ii][6] == -1)
1236  {
1237  // prism
1238  for (int i=0; i<2; i++)
1239  {
1240  K[0] = i;
1241  for (int j=0; j<2; j++)
1242  {
1243  K[1] = j;
1244  for (int k=0; k<3; k++)
1245  {
1246  K[2] = 0.5*static_cast<Real>(k);
1247  K[3] = static_cast<Real>(k % 2);
1248  basisA(Q, 6, K, H[ii], me);
1249 
1250  std::vector<Real> a1(3), a2(3), a3(3);
1251  for (int kk=0; kk<3; kk++)
1252  for (int ll=0; ll<6; ll++)
1253  {
1254  a1[kk] += Q[kk][ll]*R[cells[ii][ll]][0];
1255  a2[kk] += Q[kk][ll]*R[cells[ii][ll]][1];
1256  a3[kk] += Q[kk][ll]*R[cells[ii][ll]][2];
1257  }
1258 
1259  Real det = jac3(a1[0], a1[1], a1[2],
1260  a2[0], a2[1], a2[2],
1261  a3[0], a3[1], a3[2]);
1262  Real tr = 0;
1263  for (int kk=0; kk<3; kk++)
1264  tr += (a1[kk]*a1[kk] + a2[kk]*a2[kk] + a3[kk]*a3[kk])/3.;
1265 
1266  Real chi = 0.5*(det+std::sqrt(det*det+epsilon*epsilon));
1267  E += ((1-w)*pow(tr,1.5)/chi + 0.5*w*(v + det*det/v)/chi)/12.;
1268  if (vmin > det)
1269  vmin = det;
1270  }
1271  }
1272  }
1273 
1274  if (E > gemax)
1275  gemax = E;
1276  }
1277  else if (cells[ii][8] == -1)
1278  {
1279  // hex
1280  for (int i=0; i<2; i++)
1281  {
1282  K[0] = i;
1283  for (int j=0; j<2; j++)
1284  {
1285  K[1] = j;
1286  for (int k=0; k<2; k++)
1287  {
1288  K[2] = k;
1289  for (int l=0; l<2; l++)
1290  {
1291  K[3] = l;
1292  for (int m=0; m<2; m++)
1293  {
1294  K[4] = m;
1295  for (int nn=0; nn<2; nn++)
1296  {
1297  K[5] = nn;
1298  basisA(Q, 8, K, H[ii], me);
1299 
1300  std::vector<Real> a1(3), a2(3), a3(3);
1301  for (int kk=0; kk<3; kk++)
1302  for (int ll=0; ll<8; ll++)
1303  {
1304  a1[kk] += Q[kk][ll]*R[cells[ii][ll]][0];
1305  a2[kk] += Q[kk][ll]*R[cells[ii][ll]][1];
1306  a3[kk] += Q[kk][ll]*R[cells[ii][ll]][2];
1307  }
1308 
1309  Real det = jac3(a1[0], a1[1], a1[2],
1310  a2[0], a2[1], a2[2],
1311  a3[0], a3[1], a3[2]);
1312  Real sigma = 0.;
1313 
1314  if ((i==nn) && (j==l) && (k==m))
1315  sigma = 1./27.;
1316 
1317  if (((i==nn) && (j==l) && (k!=m)) ||
1318  ((i==nn) && (j!=l) && (k==m)) ||
1319  ((i!=nn) && (j==l) && (k==m)))
1320  sigma = 1./54.;
1321 
1322  if (((i==nn) && (j!=l) && (k!=m)) ||
1323  ((i!=nn) && (j!=l) && (k==m)) ||
1324  ((i!=nn) && (j==l) && (k!=m)))
1325  sigma = 1./108.;
1326 
1327  if ((i!=nn) && (j!=l) && (k!=m))
1328  sigma = 1./216.;
1329 
1330  Real tr = 0;
1331  for (int kk=0; kk<3; kk++)
1332  tr += (a1[kk]*a1[kk] + a2[kk]*a2[kk] + a3[kk]*a3[kk])/3.;
1333 
1334  Real chi = 0.5*(det+std::sqrt(det*det + epsilon*epsilon));
1335  E += ((1-w)*pow(tr,1.5)/chi + 0.5*w*(v + det*det/v)/chi)*sigma;
1336  if (vmin > det)
1337  vmin = det;
1338  }
1339  }
1340  }
1341  }
1342  }
1343  }
1344  if (E > gemax)
1345  gemax = E;
1346  }
1347  else
1348  {
1349  // quad tetr
1350  for (int i=0; i<4; i++)
1351  {
1352  for (int j=0; j<4; j++)
1353  {
1354  for (int k=0; k<4; k++)
1355  {
1356  switch (i)
1357  {
1358  case 0:
1359  K[0] = 0;
1360  K[1] = 0;
1361  K[2] = 0;
1362  break;
1363  case 1:
1364  K[0] = 1;
1365  K[1] = 0;
1366  K[2] = 0;
1367  break;
1368  case 2:
1369  K[0] = 0.5;
1370  K[1] = 1;
1371  K[2] = 0;
1372  break;
1373  case 3:
1374  K[0] = 0.5;
1375  K[1] = 1./3.;
1376  K[2] = 1;
1377  break;
1378  default:
1379  break;
1380  }
1381 
1382  switch (j)
1383  {
1384  case 0:
1385  K[3] = 0;
1386  K[4] = 0;
1387  K[5] = 0;
1388  break;
1389  case 1:
1390  K[3] = 1;
1391  K[4] = 0;
1392  K[5] = 0;
1393  break;
1394  case 2:
1395  K[3] = 0.5;
1396  K[4] = 1;
1397  K[5] = 0;
1398  break;
1399  case 3:
1400  K[3] = 0.5;
1401  K[4] = 1./3.;
1402  K[5] = 1;
1403  break;
1404  default:
1405  break;
1406  }
1407 
1408  switch (k)
1409  {
1410  case 0:
1411  K[6] = 0;
1412  K[7] = 0;
1413  K[8] = 0;
1414  break;
1415  case 1:
1416  K[6] = 1;
1417  K[7] = 0;
1418  K[8] = 0;
1419  break;
1420  case 2:
1421  K[6] = 0.5;
1422  K[7] = 1;
1423  K[8] = 0;
1424  break;
1425  case 3:
1426  K[6] = 0.5;
1427  K[7] = 1./3.;
1428  K[8] = 1;
1429  break;
1430  default:
1431  break;
1432  }
1433 
1434  basisA(Q, 10, K, H[ii], me);
1435 
1436  std::vector<Real> a1(3), a2(3), a3(3);
1437  for (int kk=0; kk<3; kk++)
1438  for (int ll=0; ll<10; ll++)
1439  {
1440  a1[kk] += Q[kk][ll]*R[cells[ii][ll]][0];
1441  a2[kk] += Q[kk][ll]*R[cells[ii][ll]][1];
1442  a3[kk] += Q[kk][ll]*R[cells[ii][ll]][2];
1443  }
1444 
1445  Real det = jac3(a1[0], a1[1], a1[2],
1446  a2[0], a2[1], a2[2],
1447  a3[0], a3[1], a3[2]);
1448  Real sigma = 0.;
1449 
1450  if ((i==j) && (j==k))
1451  sigma = 1./120.;
1452  else if ((i==j) || (j==k) || (i==k))
1453  sigma = 1./360.;
1454  else
1455  sigma = 1./720.;
1456 
1457  Real tr = 0;
1458  for (int kk=0; kk<3; kk++)
1459  tr += (a1[kk]*a1[kk] + a2[kk]*a2[kk] + a3[kk]*a3[kk])/3.;
1460 
1461  Real chi = 0.5*(det+std::sqrt(det*det+epsilon*epsilon));
1462  E += ((1-w)*pow(tr,1.5)/chi + 0.5*w*(v+det*det/v)/chi)*sigma;
1463  if (vmin > det)
1464  vmin = det;
1465  }
1466  }
1467  }
1468 
1469  if (E > gemax)
1470  gemax = E;
1471  }
1472  }
1473  Gamma[ii] = E;
1474  }
1475 
1476  qmin = vmin;
1477 
1478  return gemax;
1479 }

References _dim, _n_cells, basisA(), jac2(), jac3(), std::pow(), libMesh::Real, and std::sqrt().

Referenced by full_smooth().

◆ metr_data_gen()

void libMesh::VariationalMeshSmoother::metr_data_gen ( std::string  grid,
std::string  metr,
int  me 
)
private

Definition at line 3988 of file mesh_smoother_vsmoother.C.

3991 {
3992  Real det, g1, g2, g3, det_o, g1_o, g2_o, g3_o, eps=1e-3;
3993 
3994  std::vector<Real> K(9);
3995  Array2D<Real> Q(3, 3*_dim + _dim%2);
3996 
3997  // Use _mesh to determine N and ncells
3998  this->_n_nodes = _mesh.n_nodes();
3999  this->_n_cells = _mesh.n_active_elem();
4000 
4001  std::vector<int>
4002  mask(_n_nodes),
4003  mcells(_n_cells);
4004 
4005  Array2D<int> cells(_n_cells, 3*_dim + _dim%2);
4006  Array2D<Real> R(_n_nodes,_dim);
4007 
4008  readgr(R, mask, cells, mcells, mcells, mcells);
4009 
4010  // generate metric file
4011  std::ofstream metric_file(metr.c_str());
4012  if (!metric_file.good())
4013  libmesh_error_msg("Error opening metric output file.");
4014 
4015  // Use scientific notation with 6 digits
4016  metric_file << std::scientific << std::setprecision(6);
4017 
4018  int Ncells = 0;
4019  det_o = 1.;
4020  g1_o = 1.;
4021  g2_o = 1.;
4022  g3_o = 1.;
4023  for (unsigned i=0; i<_n_cells; i++)
4024  if (mcells[i] >= 0)
4025  {
4026  int nvert=0;
4027  while (cells[i][nvert] >= 0)
4028  nvert++;
4029 
4030  if (_dim == 2)
4031  {
4032  // 2D - tri and quad
4033  if (nvert == 3)
4034  {
4035  // tri
4036  basisA(Q, 3, K, Q, 1);
4037 
4038  Point a1, a2;
4039  for (int k=0; k<2; k++)
4040  for (int l=0; l<3; l++)
4041  {
4042  a1(k) += Q[k][l]*R[cells[i][l]][0];
4043  a2(k) += Q[k][l]*R[cells[i][l]][1];
4044  }
4045 
4046  det = jac2(a1(0), a1(1), a2(0), a2(1));
4047  g1 = std::sqrt(a1(0)*a1(0) + a2(0)*a2(0));
4048  g2 = std::sqrt(a1(1)*a1(1) + a2(1)*a2(1));
4049 
4050  // need to keep data from previous cell
4051  if ((std::abs(det) < eps*eps*det_o) || (det < 0))
4052  det = det_o;
4053 
4054  if ((std::abs(g1) < eps*g1_o) || (g1<0))
4055  g1 = g1_o;
4056 
4057  if ((std::abs(g2) < eps*g2_o) || (g2<0))
4058  g2 = g2_o;
4059 
4060  // write to file
4061  if (me == 2)
4062  metric_file << 1./std::sqrt(det)
4063  << " 0.000000e+00 \n0.000000e+00 "
4064  << 1./std::sqrt(det)
4065  << std::endl;
4066 
4067  if (me == 3)
4068  metric_file << 1./g1
4069  << " 0.000000e+00 \n0.000000e+00 "
4070  << 1./g2
4071  << std::endl;
4072 
4073  det_o = det;
4074  g1_o = g1;
4075  g2_o = g2;
4076  Ncells++;
4077  }
4078 
4079  if (nvert == 4)
4080  {
4081  // quad
4082 
4083  // Set up the node edge neighbor indices
4084  const unsigned node_indices[4] = {0, 1, 3, 2};
4085  const unsigned first_neighbor_indices[4] = {1, 3, 2, 0};
4086  const unsigned second_neighbor_indices[4] = {2, 0, 1, 3};
4087 
4088  // Loop over each node, compute some quantities associated
4089  // with its edge neighbors, and write them to file.
4090  for (unsigned ni=0; ni<4; ++ni)
4091  {
4092  unsigned
4093  node_index = node_indices[ni],
4094  first_neighbor_index = first_neighbor_indices[ni],
4095  second_neighbor_index = second_neighbor_indices[ni];
4096 
4097  Real
4098  node_x = R[cells[i][node_index]][0],
4099  node_y = R[cells[i][node_index]][1],
4100  first_neighbor_x = R[cells[i][first_neighbor_index]][0],
4101  first_neighbor_y = R[cells[i][first_neighbor_index]][1],
4102  second_neighbor_x = R[cells[i][second_neighbor_index]][0],
4103  second_neighbor_y = R[cells[i][second_neighbor_index]][1];
4104 
4105 
4106  // "dx"
4107  Point a1(first_neighbor_x - node_x,
4108  second_neighbor_x - node_x);
4109 
4110  // "dy"
4111  Point a2(first_neighbor_y - node_y,
4112  second_neighbor_y - node_y);
4113 
4114  det = jac2(a1(0), a1(1), a2(0), a2(1));
4115  g1 = std::sqrt(a1(0)*a1(0) + a2(0)*a2(0));
4116  g2 = std::sqrt(a1(1)*a1(1) + a2(1)*a2(1));
4117 
4118  // need to keep data from previous cell
4119  if ((std::abs(det) < eps*eps*det_o) || (det < 0))
4120  det = det_o;
4121 
4122  if ((std::abs(g1) < eps*g1_o) || (g1 < 0))
4123  g1 = g1_o;
4124 
4125  if ((std::abs(g2) < eps*g2_o) || (g2 < 0))
4126  g2 = g2_o;
4127 
4128  // write to file
4129  if (me == 2)
4130  metric_file << 1./std::sqrt(det) << " "
4131  << 0.5/std::sqrt(det) << " \n0.000000e+00 "
4132  << 0.5*std::sqrt(3./det)
4133  << std::endl;
4134 
4135  if (me == 3)
4136  metric_file << 1./g1 << " "
4137  << 0.5/g2 << "\n0.000000e+00 "
4138  << 0.5*std::sqrt(3.)/g2
4139  << std::endl;
4140 
4141  det_o = det;
4142  g1_o = g1;
4143  g2_o = g2;
4144  Ncells++;
4145  }
4146  } // end QUAD case
4147  } // end _dim==2
4148 
4149 #if LIBMESH_DIM > 2
4150  if (_dim == 3)
4151  {
4152  // 3D - tetr and hex
4153 
4154  if (nvert == 4)
4155  {
4156  // tetr
4157  basisA(Q, 4, K, Q, 1);
4158 
4159  Point a1, a2, a3;
4160 
4161  for (int k=0; k<3; k++)
4162  for (int l=0; l<4; l++)
4163  {
4164  a1(k) += Q[k][l]*R[cells[i][l]][0];
4165  a2(k) += Q[k][l]*R[cells[i][l]][1];
4166  a3(k) += Q[k][l]*R[cells[i][l]][2];
4167  }
4168 
4169  det = jac3(a1(0), a1(1), a1(2),
4170  a2(0), a2(1), a2(2),
4171  a3(0), a3(1), a3(2));
4172  g1 = std::sqrt(a1(0)*a1(0) + a2(0)*a2(0) + a3(0)*a3(0));
4173  g2 = std::sqrt(a1(1)*a1(1) + a2(1)*a2(1) + a3(1)*a3(1));
4174  g3 = std::sqrt(a1(2)*a1(2) + a2(2)*a2(2) + a3(2)*a3(2));
4175 
4176  // need to keep data from previous cell
4177  if ((std::abs(det) < eps*eps*eps*det_o) || (det < 0))
4178  det = det_o;
4179 
4180  if ((std::abs(g1) < eps*g1_o) || (g1 < 0))
4181  g1 = g1_o;
4182 
4183  if ((std::abs(g2) < eps*g2_o) || (g2 < 0))
4184  g2 = g2_o;
4185 
4186  if ((std::abs(g3) < eps*g3_o) || (g3 < 0))
4187  g3 = g3_o;
4188 
4189  // write to file
4190  if (me == 2)
4191  metric_file << 1./pow(det, 1./3.)
4192  << " 0.000000e+00 0.000000e+00\n0.000000e+00 "
4193  << 1./pow(det, 1./3.)
4194  << " 0.000000e+00\n0.000000e+00 0.000000e+00 "
4195  << 1./pow(det, 1./3.)
4196  << std::endl;
4197 
4198  if (me == 3)
4199  metric_file << 1./g1
4200  << " 0.000000e+00 0.000000e+00\n0.000000e+00 "
4201  << 1./g2
4202  << " 0.000000e+00\n0.000000e+00 0.000000e+00 "
4203  << 1./g3
4204  << std::endl;
4205 
4206  det_o = det;
4207  g1_o = g1;
4208  g2_o = g2;
4209  g3_o = g3;
4210  Ncells++;
4211  }
4212 
4213  if (nvert == 8)
4214  {
4215  // hex
4216 
4217  // Set up the node edge neighbor indices
4218  const unsigned node_indices[8] = {0, 1, 3, 2, 4, 5, 7, 6};
4219  const unsigned first_neighbor_indices[8] = {1, 3, 2, 0, 6, 4, 5, 7};
4220  const unsigned second_neighbor_indices[8] = {2, 0, 1, 3, 5, 7, 6, 4};
4221  const unsigned third_neighbor_indices[8] = {4, 5, 7, 6, 0, 1, 3, 2};
4222 
4223  // Loop over each node, compute some quantities associated
4224  // with its edge neighbors, and write them to file.
4225  for (unsigned ni=0; ni<8; ++ni)
4226  {
4227  unsigned
4228  node_index = node_indices[ni],
4229  first_neighbor_index = first_neighbor_indices[ni],
4230  second_neighbor_index = second_neighbor_indices[ni],
4231  third_neighbor_index = third_neighbor_indices[ni];
4232 
4233  Real
4234  node_x = R[cells[i][node_index]][0],
4235  node_y = R[cells[i][node_index]][1],
4236  node_z = R[cells[i][node_index]][2],
4237  first_neighbor_x = R[cells[i][first_neighbor_index]][0],
4238  first_neighbor_y = R[cells[i][first_neighbor_index]][1],
4239  first_neighbor_z = R[cells[i][first_neighbor_index]][2],
4240  second_neighbor_x = R[cells[i][second_neighbor_index]][0],
4241  second_neighbor_y = R[cells[i][second_neighbor_index]][1],
4242  second_neighbor_z = R[cells[i][second_neighbor_index]][2],
4243  third_neighbor_x = R[cells[i][third_neighbor_index]][0],
4244  third_neighbor_y = R[cells[i][third_neighbor_index]][1],
4245  third_neighbor_z = R[cells[i][third_neighbor_index]][2];
4246 
4247  // "dx"
4248  Point a1(first_neighbor_x - node_x,
4249  second_neighbor_x - node_x,
4250  third_neighbor_x - node_x);
4251 
4252  // "dy"
4253  Point a2(first_neighbor_y - node_y,
4254  second_neighbor_y - node_y,
4255  third_neighbor_y - node_y);
4256 
4257  // "dz"
4258  Point a3(first_neighbor_z - node_z,
4259  second_neighbor_z - node_z,
4260  third_neighbor_z - node_z);
4261 
4262  det = jac3(a1(0), a1(1), a1(2),
4263  a2(0), a2(1), a2(2),
4264  a3(0), a3(1), a3(2));
4265  g1 = std::sqrt(a1(0)*a1(0) + a2(0)*a2(0) + a3(0)*a3(0));
4266  g2 = std::sqrt(a1(1)*a1(1) + a2(1)*a2(1) + a3(1)*a3(1));
4267  g3 = std::sqrt(a1(2)*a1(2) + a2(2)*a2(2) + a3(2)*a3(2));
4268 
4269  // need to keep data from previous cell
4270  if ((std::abs(det) < eps*eps*eps*det_o) || (det < 0))
4271  det = det_o;
4272 
4273  if ((std::abs(g1) < eps*g1_o) || (g1 < 0))
4274  g1 = g1_o;
4275 
4276  if ((std::abs(g2) < eps*g2_o) || (g2 < 0))
4277  g2 = g2_o;
4278 
4279  if ((std::abs(g3) < eps*g3_o) || (g3 < 0))
4280  g3 = g3_o;
4281 
4282  // write to file
4283  if (me == 2)
4284  metric_file << 1./pow(det, 1./3.) << " "
4285  << 0.5/pow(det, 1./3.) << " "
4286  << 0.5/pow(det, 1./3.) << "\n0.000000e+00 "
4287  << 0.5*std::sqrt(3.)/pow(det, 1./3.) << " "
4288  << 0.5/(std::sqrt(3.)*pow(det, 1./3.)) << "\n0.000000e+00 0.000000e+00 "
4289  << std::sqrt(2/3.)/pow(det, 1./3.)
4290  << std::endl;
4291 
4292  if (me == 3)
4293  metric_file << 1./g1 << " "
4294  << 0.5/g2 << " "
4295  << 0.5/g3 << "\n0.000000e+00 "
4296  << 0.5*std::sqrt(3.)/g2 << " "
4297  << 0.5/(std::sqrt(3.)*g3) << "\n0.000000e+00 0.000000e+00 "
4298  << std::sqrt(2./3.)/g3
4299  << std::endl;
4300 
4301  det_o = det;
4302  g1_o = g1;
4303  g2_o = g2;
4304  g3_o = g3;
4305  Ncells++;
4306  } // end for ni
4307  } // end hex
4308  } // end dim==3
4309 #endif // LIBMESH_DIM > 2
4310  }
4311 
4312  // Done with the metric file
4313  metric_file.close();
4314 
4315  std::ofstream grid_file(grid.c_str());
4316  if (!grid_file.good())
4317  libmesh_error_msg("Error opening file: " << grid);
4318 
4319  grid_file << _dim << "\n" << _n_nodes << "\n" << Ncells << "\n0" << std::endl;
4320 
4321  // Use scientific notation with 6 digits
4322  grid_file << std::scientific << std::setprecision(6);
4323 
4324  for (dof_id_type i=0; i<_n_nodes; i++)
4325  {
4326  // node coordinates
4327  for (unsigned j=0; j<_dim; j++)
4328  grid_file << R[i][j] << " ";
4329 
4330  grid_file << mask[i] << std::endl;
4331  }
4332 
4333  for (dof_id_type i=0; i<_n_cells; i++)
4334  if (mcells[i] >= 0)
4335  {
4336  // cell connectivity
4337  int nvert = 0;
4338  while (cells[i][nvert] >= 0)
4339  nvert++;
4340 
4341  if ((nvert == 3) || ((_dim == 3) && (nvert == 4)))
4342  {
4343  // tri & tetr
4344  for (int j=0; j<nvert; j++)
4345  grid_file << cells[i][j] << " ";
4346 
4347  for (unsigned j=nvert; j<3*_dim + _dim%2; j++)
4348  grid_file << "-1 ";
4349 
4350  grid_file << mcells[i] << std::endl;
4351  }
4352 
4353  if ((_dim == 2) && (nvert == 4))
4354  {
4355  // quad
4356  grid_file << cells[i][0] << " " << cells[i][1] << " " << cells[i][2] << " -1 -1 -1 0" << std::endl;
4357  grid_file << cells[i][1] << " " << cells[i][3] << " " << cells[i][0] << " -1 -1 -1 0" << std::endl;
4358  grid_file << cells[i][3] << " " << cells[i][2] << " " << cells[i][1] << " -1 -1 -1 0" << std::endl;
4359  grid_file << cells[i][2] << " " << cells[i][0] << " " << cells[i][3] << " -1 -1 -1 0" << std::endl;
4360  }
4361 
4362  if (nvert == 8)
4363  {
4364  // hex
4365  grid_file << cells[i][0] << " " << cells[i][1] << " " << cells[i][2] << " " << cells[i][4] << " -1 -1 -1 -1 -1 -1 0" << std::endl;
4366  grid_file << cells[i][1] << " " << cells[i][3] << " " << cells[i][0] << " " << cells[i][5] << " -1 -1 -1 -1 -1 -1 0" << std::endl;
4367  grid_file << cells[i][3] << " " << cells[i][2] << " " << cells[i][1] << " " << cells[i][7] << " -1 -1 -1 -1 -1 -1 0" << std::endl;
4368  grid_file << cells[i][2] << " " << cells[i][0] << " " << cells[i][3] << " " << cells[i][6] << " -1 -1 -1 -1 -1 -1 0" << std::endl;
4369  grid_file << cells[i][4] << " " << cells[i][6] << " " << cells[i][5] << " " << cells[i][0] << " -1 -1 -1 -1 -1 -1 0" << std::endl;
4370  grid_file << cells[i][5] << " " << cells[i][4] << " " << cells[i][7] << " " << cells[i][1] << " -1 -1 -1 -1 -1 -1 0" << std::endl;
4371  grid_file << cells[i][7] << " " << cells[i][5] << " " << cells[i][6] << " " << cells[i][3] << " -1 -1 -1 -1 -1 -1 0" << std::endl;
4372  grid_file << cells[i][6] << " " << cells[i][7] << " " << cells[i][4] << " " << cells[i][2] << " -1 -1 -1 -1 -1 -1 0" << std::endl;
4373  }
4374  }
4375 }

References _dim, libMesh::MeshSmoother::_mesh, _n_cells, _n_nodes, std::abs(), basisA(), jac2(), jac3(), libMesh::MeshBase::n_active_elem(), libMesh::MeshBase::n_nodes(), std::pow(), readgr(), libMesh::Real, and std::sqrt().

Referenced by smooth().

◆ minJ()

Real libMesh::VariationalMeshSmoother::minJ ( Array2D< Real > &  R,
const std::vector< int > &  mask,
const Array2D< int > &  cells,
const std::vector< int > &  mcells,
Real  epsilon,
Real  w,
int  me,
const Array3D< Real > &  H,
Real  vol,
const std::vector< int > &  edges,
const std::vector< int > &  hnodes,
int  msglev,
Real &  Vmin,
Real &  emax,
Real &  qmin,
int  adp,
const std::vector< Real > &  afun 
)
private

Definition at line 1887 of file mesh_smoother_vsmoother.C.

1904 {
1905  // columns - max number of nonzero entries in every row of global matrix
1906  int columns = _dim*_dim*10;
1907 
1908  // local Hessian matrix
1909  Array3D<Real> W(_dim, 3*_dim + _dim%2, 3*_dim + _dim%2);
1910 
1911  // F - local gradient
1912  Array2D<Real> F(_dim, 3*_dim + _dim%2);
1913 
1914  Array2D<Real> Rpr(_n_nodes, _dim);
1915 
1916  // P - minimization direction
1917  Array2D<Real> P(_n_nodes, _dim);
1918 
1919  // A, JA - internal form of global matrix
1920  Array2D<int> JA(_dim*_n_nodes, columns);
1921  Array2D<Real> A(_dim*_n_nodes, columns);
1922 
1923  // G - adaptation metric
1924  Array2D<Real> G(_n_cells, _dim);
1925 
1926  // rhs for solver
1927  std::vector<Real> b(_dim*_n_nodes);
1928 
1929  // u - solution vector
1930  std::vector<Real> u(_dim*_n_nodes);
1931 
1932  // matrix
1933  std::vector<Real> a(_dim*_n_nodes*columns);
1934  std::vector<int> ia(_dim*_n_nodes + 1);
1935  std::vector<int> ja(_dim*_n_nodes*columns);
1936 
1937  // nonzero - norm of gradient
1938  Real nonzero = 0.;
1939 
1940  // Jpr - value of functional
1941  Real Jpr = 0.;
1942 
1943  // find minimization direction P
1944  for (dof_id_type i=0; i<_n_cells; i++)
1945  {
1946  int nvert = 0;
1947  while (cells[i][nvert] >= 0)
1948  nvert++;
1949 
1950  // determination of local matrices on each cell
1951  for (unsigned j=0; j<_dim; j++)
1952  {
1953  G[i][j] = 0; // adaptation metric G is held constant throughout minJ run
1954  if (adp < 0)
1955  {
1956  for (auto k : IntRange<int>(0, std::abs(adp)))
1957  G[i][j] += afun[i*(-adp)+k]; // cell-based adaptivity is computed here
1958  }
1959  }
1960  for (unsigned index=0; index<_dim; index++)
1961  {
1962  // initialize local matrices
1963  for (unsigned k=0; k<3*_dim + _dim%2; k++)
1964  {
1965  F[index][k] = 0;
1966 
1967  for (unsigned j=0; j<3*_dim + _dim%2; j++)
1968  W[index][k][j] = 0;
1969  }
1970  }
1971  if (mcells[i] >= 0)
1972  {
1973  // if cell is not excluded
1974  Real lVmin, lqmin;
1975  Jpr += localP(W, F, R, cells[i], mask, epsilon, w, nvert, H[i],
1976  me, vol, 0, lVmin, lqmin, adp, afun, G[i]);
1977  }
1978  else
1979  {
1980  for (unsigned index=0; index<_dim; index++)
1981  for (int j=0; j<nvert; j++)
1982  W[index][j][j] = 1;
1983  }
1984 
1985  // assembly of an upper triangular part of a global matrix A
1986  for (unsigned index=0; index<_dim; index++)
1987  {
1988  for (int l=0; l<nvert; l++)
1989  {
1990  for (int m=0; m<nvert; m++)
1991  {
1992  if ((W[index][l][m] != 0) &&
1993  (cells[i][m] >= cells[i][l]))
1994  {
1995  int sch = 0;
1996  int ind = 1;
1997  while (ind != 0)
1998  {
1999  if (A[cells[i][l] + index*_n_nodes][sch] != 0)
2000  {
2001  if (JA[cells[i][l] + index*_n_nodes][sch] == static_cast<int>(cells[i][m] + index*_n_nodes))
2002  {
2003  A[cells[i][l] + index*_n_nodes][sch] = A[cells[i][l] + index*_n_nodes][sch] + W[index][l][m];
2004  ind=0;
2005  }
2006  else
2007  sch++;
2008  }
2009  else
2010  {
2011  A[cells[i][l] + index*_n_nodes][sch] = W[index][l][m];
2012  JA[cells[i][l] + index*_n_nodes][sch] = cells[i][m] + index*_n_nodes;
2013  ind = 0;
2014  }
2015 
2016  if (sch > columns-1)
2017  _logfile << "error: # of nonzero entries in the "
2018  << cells[i][l]
2019  << " row of Hessian ="
2020  << sch
2021  << ">= columns="
2022  << columns
2023  << std::endl;
2024  }
2025  }
2026  }
2027  b[cells[i][l] + index*_n_nodes] = b[cells[i][l] + index*_n_nodes] - F[index][l];
2028  }
2029  }
2030  // end of matrix A
2031  }
2032 
2033  // HN correction
2034 
2035  // tolerance for HN being the mid-edge node for its parents
2036  Real Tau_hn = pow(vol, 1./static_cast<Real>(_dim))*1e-10;
2037  for (dof_id_type i=0; i<_n_hanging_edges; i++)
2038  {
2039  int ind_i = hnodes[i];
2040  int ind_j = edges[2*i];
2041  int ind_k = edges[2*i+1];
2042 
2043  for (unsigned j=0; j<_dim; j++)
2044  {
2045  int g_i = int(R[ind_i][j] - 0.5*(R[ind_j][j]+R[ind_k][j]));
2046  Jpr += g_i*g_i/(2*Tau_hn);
2047  b[ind_i + j*_n_nodes] -= g_i/Tau_hn;
2048  b[ind_j + j*_n_nodes] += 0.5*g_i/Tau_hn;
2049  b[ind_k + j*_n_nodes] += 0.5*g_i/Tau_hn;
2050  }
2051 
2052  for (int ind=0; ind<columns; ind++)
2053  {
2054  if (JA[ind_i][ind] == ind_i)
2055  A[ind_i][ind] += 1./Tau_hn;
2056 
2057  if (JA[ind_i+_n_nodes][ind] == static_cast<int>(ind_i+_n_nodes))
2058  A[ind_i+_n_nodes][ind] += 1./Tau_hn;
2059 
2060  if (_dim == 3)
2061  if (JA[ind_i+2*_n_nodes][ind] == static_cast<int>(ind_i+2*_n_nodes))
2062  A[ind_i+2*_n_nodes][ind] += 1./Tau_hn;
2063 
2064  if ((JA[ind_i][ind] == ind_j) ||
2065  (JA[ind_i][ind] == ind_k))
2066  A[ind_i][ind] -= 0.5/Tau_hn;
2067 
2068  if ((JA[ind_i+_n_nodes][ind] == static_cast<int>(ind_j+_n_nodes)) ||
2069  (JA[ind_i+_n_nodes][ind] == static_cast<int>(ind_k+_n_nodes)))
2070  A[ind_i+_n_nodes][ind] -= 0.5/Tau_hn;
2071 
2072  if (_dim == 3)
2073  if ((JA[ind_i+2*_n_nodes][ind] == static_cast<int>(ind_j+2*_n_nodes)) ||
2074  (JA[ind_i+2*_n_nodes][ind] == static_cast<int>(ind_k+2*_n_nodes)))
2075  A[ind_i+2*_n_nodes][ind] -= 0.5/Tau_hn;
2076 
2077  if (JA[ind_j][ind] == ind_i)
2078  A[ind_j][ind] -= 0.5/Tau_hn;
2079 
2080  if (JA[ind_k][ind] == ind_i)
2081  A[ind_k][ind] -= 0.5/Tau_hn;
2082 
2083  if (JA[ind_j+_n_nodes][ind] == static_cast<int>(ind_i+_n_nodes))
2084  A[ind_j+_n_nodes][ind] -= 0.5/Tau_hn;
2085 
2086  if (JA[ind_k+_n_nodes][ind] == static_cast<int>(ind_i+_n_nodes))
2087  A[ind_k+_n_nodes][ind] -= 0.5/Tau_hn;
2088 
2089  if (_dim == 3)
2090  if (JA[ind_j+2*_n_nodes][ind] == static_cast<int>(ind_i+2*_n_nodes))
2091  A[ind_j+2*_n_nodes][ind] -= 0.5/Tau_hn;
2092 
2093  if (_dim == 3)
2094  if (JA[ind_k+2*_n_nodes][ind] == static_cast<int>(ind_i+2*_n_nodes))
2095  A[ind_k+2*_n_nodes][ind] -= 0.5/Tau_hn;
2096 
2097  if ((JA[ind_j][ind] == ind_j) ||
2098  (JA[ind_j][ind] == ind_k))
2099  A[ind_j][ind] += 0.25/Tau_hn;
2100 
2101  if ((JA[ind_k][ind] == ind_j) ||
2102  (JA[ind_k][ind] == ind_k))
2103  A[ind_k][ind] += 0.25/Tau_hn;
2104 
2105  if ((JA[ind_j+_n_nodes][ind] == static_cast<int>(ind_j+_n_nodes)) ||
2106  (JA[ind_j+_n_nodes][ind] == static_cast<int>(ind_k+_n_nodes)))
2107  A[ind_j+_n_nodes][ind] += 0.25/Tau_hn;
2108 
2109  if ((JA[ind_k+_n_nodes][ind] == static_cast<int>(ind_j+_n_nodes)) ||
2110  (JA[ind_k+_n_nodes][ind] == static_cast<int>(ind_k+_n_nodes)))
2111  A[ind_k+_n_nodes][ind] += 0.25/Tau_hn;
2112 
2113  if (_dim == 3)
2114  if ((JA[ind_j+2*_n_nodes][ind] == static_cast<int>(ind_j+2*_n_nodes)) ||
2115  (JA[ind_j+2*_n_nodes][ind] == static_cast<int>(ind_k+2*_n_nodes)))
2116  A[ind_j+2*_n_nodes][ind] += 0.25/Tau_hn;
2117 
2118  if (_dim==3)
2119  if ((JA[ind_k+2*_n_nodes][ind] == static_cast<int>(ind_j+2*_n_nodes)) ||
2120  (JA[ind_k+2*_n_nodes][ind] == static_cast<int>(ind_k+2*_n_nodes)))
2121  A[ind_k+2*_n_nodes][ind] += 0.25/Tau_hn;
2122  }
2123  }
2124 
2125  // ||\grad J||_2
2126  for (dof_id_type i=0; i<_dim*_n_nodes; i++)
2127  nonzero += b[i]*b[i];
2128 
2129  // sort matrix A
2130  for (dof_id_type i=0; i<_dim*_n_nodes; i++)
2131  {
2132  for (int j=columns-1; j>1; j--)
2133  {
2134  for (int k=0; k<j; k++)
2135  {
2136  if (JA[i][k] > JA[i][k+1])
2137  {
2138  int sch = JA[i][k];
2139  JA[i][k] = JA[i][k+1];
2140  JA[i][k+1] = sch;
2141  Real tmp = A[i][k];
2142  A[i][k] = A[i][k+1];
2143  A[i][k+1] = tmp;
2144  }
2145  }
2146  }
2147  }
2148 
2149  Real eps = std::sqrt(vol)*1e-9;
2150 
2151  // solver for P (unconstrained)
2152  ia[0] = 0;
2153  {
2154  int j = 0;
2155  for (dof_id_type i=0; i<_dim*_n_nodes; i++)
2156  {
2157  u[i] = 0;
2158  int nz = 0;
2159  for (int k=0; k<columns; k++)
2160  {
2161  if (A[i][k] != 0)
2162  {
2163  nz++;
2164  ja[j] = JA[i][k]+1;
2165  a[j] = A[i][k];
2166  j++;
2167  }
2168  }
2169  ia[i+1] = ia[i] + nz;
2170  }
2171  }
2172 
2173  dof_id_type m = _dim*_n_nodes;
2174  int sch = (msglev >= 3) ? 1 : 0;
2175 
2176  solver(m, ia, ja, a, u, b, eps, 100, sch);
2177  // sol_pcg_pj(m, ia, ja, a, u, b, eps, 100, sch);
2178 
2179  for (dof_id_type i=0; i<_n_nodes; i++)
2180  {
2181  //ensure fixed nodes are not moved
2182  for (unsigned index=0; index<_dim; index++)
2183  if (mask[i] == 1)
2184  P[i][index] = 0;
2185  else
2186  P[i][index] = u[i+index*_n_nodes];
2187  }
2188 
2189  // P is determined
2190  if (msglev >= 4)
2191  {
2192  for (dof_id_type i=0; i<_n_nodes; i++)
2193  {
2194  for (unsigned j=0; j<_dim; j++)
2195  if (P[i][j] != 0)
2196  _logfile << "P[" << i << "][" << j << "]=" << P[i][j];
2197  }
2198  }
2199 
2200  // local minimization problem, determination of tau
2201  if (msglev >= 3)
2202  _logfile << "dJ=" << std::sqrt(nonzero) << " J0=" << Jpr << std::endl;
2203 
2204  Real
2205  J = 1.e32,
2206  tau = 0.,
2207  gVmin = 0.,
2208  gemax = 0.,
2209  gqmin = 0.;
2210 
2211  int j = 1;
2212 
2213  while ((Jpr <= J) && (j > -30))
2214  {
2215  j = j-1;
2216 
2217  // search through finite # of values for tau
2218  tau = pow(2., j);
2219  for (dof_id_type i=0; i<_n_nodes; i++)
2220  for (unsigned k=0; k<_dim; k++)
2221  Rpr[i][k] = R[i][k] + tau*P[i][k];
2222 
2223  J = 0;
2224  gVmin = 1e32;
2225  gemax = -1e32;
2226  gqmin = 1e32;
2227  for (dof_id_type i=0; i<_n_cells; i++)
2228  {
2229  if (mcells[i] >= 0)
2230  {
2231  int nvert = 0;
2232  while (cells[i][nvert] >= 0)
2233  nvert++;
2234 
2235  Real lVmin, lqmin;
2236  Real lemax = localP(W, F, Rpr, cells[i], mask, epsilon, w, nvert, H[i], me, vol, 1, lVmin, lqmin, adp, afun, G[i]);
2237 
2238  J += lemax;
2239  if (gVmin > lVmin)
2240  gVmin = lVmin;
2241 
2242  if (gemax < lemax)
2243  gemax = lemax;
2244 
2245  if (gqmin > lqmin)
2246  gqmin = lqmin;
2247 
2248  // HN correction
2249  for (dof_id_type ii=0; ii<_n_hanging_edges; ii++)
2250  {
2251  int ind_i = hnodes[ii];
2252  int ind_j = edges[2*ii];
2253  int ind_k = edges[2*ii+1];
2254  for (unsigned jj=0; jj<_dim; jj++)
2255  {
2256  int g_i = int(Rpr[ind_i][jj] - 0.5*(Rpr[ind_j][jj]+Rpr[ind_k][jj]));
2257  J += g_i*g_i/(2*Tau_hn);
2258  }
2259  }
2260  }
2261  }
2262  if (msglev >= 3)
2263  _logfile << "tau=" << tau << " J=" << J << std::endl;
2264  }
2265 
2266 
2267  Real
2268  T = 0.,
2269  gtmin0 = 0.,
2270  gtmax0 = 0.,
2271  gqmin0 = 0.;
2272 
2273  if (j != -30)
2274  {
2275  Jpr = J;
2276  for (unsigned i=0; i<_n_nodes; i++)
2277  for (unsigned k=0; k<_dim; k++)
2278  Rpr[i][k] = R[i][k] + tau*0.5*P[i][k];
2279 
2280  J = 0;
2281  gtmin0 = 1e32;
2282  gtmax0 = -1e32;
2283  gqmin0 = 1e32;
2284  for (dof_id_type i=0; i<_n_cells; i++)
2285  {
2286  if (mcells[i] >= 0)
2287  {
2288  int nvert = 0;
2289  while (cells[i][nvert] >= 0)
2290  nvert++;
2291 
2292  Real lVmin, lqmin;
2293  Real lemax = localP(W, F, Rpr, cells[i], mask, epsilon, w, nvert, H[i], me, vol, 1, lVmin, lqmin, adp, afun, G[i]);
2294  J += lemax;
2295 
2296  if (gtmin0 > lVmin)
2297  gtmin0 = lVmin;
2298 
2299  if (gtmax0 < lemax)
2300  gtmax0 = lemax;
2301 
2302  if (gqmin0 > lqmin)
2303  gqmin0 = lqmin;
2304 
2305  // HN correction
2306  for (dof_id_type ii=0; ii<_n_hanging_edges; ii++)
2307  {
2308  int ind_i = hnodes[ii];
2309  int ind_j = edges[2*ii];
2310  int ind_k = edges[2*ii+1];
2311 
2312  for (unsigned jj=0; jj<_dim; jj++)
2313  {
2314  int g_i = int(Rpr[ind_i][jj] - 0.5*(Rpr[ind_j][jj] + Rpr[ind_k][jj]));
2315  J += g_i*g_i/(2*Tau_hn);
2316  }
2317  }
2318  }
2319  }
2320  }
2321 
2322  if (Jpr > J)
2323  {
2324  T = 0.5*tau;
2325  // Set up return values passed by reference
2326  Vmin = gtmin0;
2327  emax = gtmax0;
2328  qmin = gqmin0;
2329  }
2330  else
2331  {
2332  T = tau;
2333  J = Jpr;
2334  // Set up return values passed by reference
2335  Vmin = gVmin;
2336  emax = gemax;
2337  qmin = gqmin;
2338  }
2339 
2340  nonzero = 0;
2341  for (dof_id_type j2=0; j2<_n_nodes; j2++)
2342  for (unsigned k=0; k<_dim; k++)
2343  {
2344  R[j2][k] = R[j2][k] + T*P[j2][k];
2345  nonzero += T*P[j2][k]*T*P[j2][k];
2346  }
2347 
2348  if (msglev >= 2)
2349  _logfile << "tau=" << T << ", J=" << J << std::endl;
2350 
2351  return std::sqrt(nonzero);
2352 }

References _dim, _logfile, _n_cells, _n_hanging_edges, _n_nodes, A, std::abs(), int, localP(), std::pow(), libMesh::Real, solver(), and std::sqrt().

Referenced by full_smooth().

◆ minJ_BC()

Real libMesh::VariationalMeshSmoother::minJ_BC ( Array2D< Real > &  R,
const std::vector< int > &  mask,
const Array2D< int > &  cells,
const std::vector< int > &  mcells,
Real  epsilon,
Real  w,
int  me,
const Array3D< Real > &  H,
Real  vol,
int  msglev,
Real &  Vmin,
Real &  emax,
Real &  qmin,
int  adp,
const std::vector< Real > &  afun,
int  NCN 
)
private

Definition at line 2359 of file mesh_smoother_vsmoother.C.

2375 {
2376  // new form of matrices, 5 iterations for minL
2377  Real tau = 0., J = 0., T, Jpr, L, lVmin, lqmin, gVmin = 0., gqmin = 0., gVmin0 = 0.,
2378  gqmin0 = 0., lemax, gemax = 0., gemax0 = 0.;
2379 
2380  // array of sliding BN
2381  std::vector<int> Bind(NCN);
2382  std::vector<Real> lam(NCN);
2383  std::vector<Real> hm(2*_n_nodes);
2384  std::vector<Real> Plam(NCN);
2385 
2386  // holds constraints = local approximation to the boundary
2387  std::vector<Real> constr(4*NCN);
2388 
2389  Array2D<Real> F(2, 6);
2390  Array3D<Real> W(2, 6, 6);
2391  Array2D<Real> Rpr(_n_nodes, 2);
2392  Array2D<Real> P(_n_nodes, 2);
2393 
2394  std::vector<Real> b(2*_n_nodes);
2395 
2396  Array2D<Real> G(_n_cells, 6);
2397 
2398  // assembler of constraints
2399  const Real eps = std::sqrt(vol)*1e-9;
2400 
2401  for (int i=0; i<4*NCN; i++)
2402  constr[i] = 1./eps;
2403 
2404  {
2405  int I = 0;
2406  for (dof_id_type i=0; i<_n_nodes; i++)
2407  if (mask[i] == 2)
2408  {
2409  Bind[I] = i;
2410  I++;
2411  }
2412  }
2413 
2414  for (int I=0; I<NCN; I++)
2415  {
2416  // The boundary connectivity loop sets up the j and k indices
2417  // but I am not sure about the logic of this: j and k are overwritten
2418  // at every iteration of the boundary connectivity loop and only used
2419  // *after* the boundary connectivity loop -- this seems like a bug but
2420  // I maintained the original behavior of the algorithm [JWP].
2421  int
2422  i = Bind[I],
2423  j = 0,
2424  k = 0,
2425  ind = 0;
2426 
2427  // boundary connectivity
2428  for (dof_id_type l=0; l<_n_cells; l++)
2429  {
2430  int nvert = 0;
2431 
2432  while (cells[l][nvert] >= 0)
2433  nvert++;
2434 
2435  switch (nvert)
2436  {
2437  case 3: // tri
2438  if (i == cells[l][0])
2439  {
2440  if (mask[cells[l][1]] > 0)
2441  {
2442  if (ind == 0)
2443  j = cells[l][1];
2444  else
2445  k = cells[l][1];
2446  ind++;
2447  }
2448 
2449  if (mask[cells[l][2]] > 0)
2450  {
2451  if (ind == 0)
2452  j = cells[l][2];
2453  else
2454  k = cells[l][2];
2455  ind++;
2456  }
2457  }
2458 
2459  if (i == cells[l][1])
2460  {
2461  if (mask[cells[l][0]] > 0)
2462  {
2463  if (ind == 0)
2464  j = cells[l][0];
2465  else
2466  k = cells[l][0];
2467  ind++;
2468  }
2469 
2470  if (mask[cells[l][2]] > 0)
2471  {
2472  if (ind == 0)
2473  j = cells[l][2];
2474  else
2475  k = cells[l][2];
2476  ind++;
2477  }
2478  }
2479 
2480  if (i == cells[l][2])
2481  {
2482  if (mask[cells[l][1]] > 0)
2483  {
2484  if (ind == 0)
2485  j = cells[l][1];
2486  else
2487  k = cells[l][1];
2488  ind++;
2489  }
2490 
2491  if (mask[cells[l][0]] > 0)
2492  {
2493  if (ind == 0)
2494  j = cells[l][0];
2495  else
2496  k = cells[l][0];
2497  ind++;
2498  }
2499  }
2500  break;
2501 
2502  case 4: // quad
2503  if ((i == cells[l][0]) ||
2504  (i == cells[l][3]))
2505  {
2506  if (mask[cells[l][1]] > 0)
2507  {
2508  if (ind == 0)
2509  j = cells[l][1];
2510  else
2511  k = cells[l][1];
2512  ind++;
2513  }
2514 
2515  if (mask[cells[l][2]] > 0)
2516  {
2517  if (ind == 0)
2518  j = cells[l][2];
2519  else
2520  k = cells[l][2];
2521  ind++;
2522  }
2523  }
2524 
2525  if ((i == cells[l][1]) ||
2526  (i == cells[l][2]))
2527  {
2528  if (mask[cells[l][0]] > 0)
2529  {
2530  if (ind == 0)
2531  j = cells[l][0];
2532  else
2533  k = cells[l][0];
2534  ind++;
2535  }
2536 
2537  if (mask[cells[l][3]] > 0)
2538  {
2539  if (ind == 0)
2540  j = cells[l][3];
2541  else
2542  k = cells[l][3];
2543  ind++;
2544  }
2545  }
2546  break;
2547 
2548  case 6: //quad tri
2549  if (i == cells[l][0])
2550  {
2551  if (mask[cells[l][1]] > 0)
2552  {
2553  if (ind == 0)
2554  j = cells[l][5];
2555  else
2556  k = cells[l][5];
2557  ind++;
2558  }
2559 
2560  if (mask[cells[l][2]] > 0)
2561  {
2562  if (ind == 0)
2563  j = cells[l][4];
2564  else
2565  k = cells[l][4];
2566  ind++;
2567  }
2568  }
2569 
2570  if (i == cells[l][1])
2571  {
2572  if (mask[cells[l][0]] > 0)
2573  {
2574  if (ind == 0)
2575  j = cells[l][5];
2576  else
2577  k = cells[l][5];
2578  ind++;
2579  }
2580 
2581  if (mask[cells[l][2]] > 0)
2582  {
2583  if (ind == 0)
2584  j = cells[l][3];
2585  else
2586  k = cells[l][3];
2587  ind++;
2588  }
2589  }
2590 
2591  if (i == cells[l][2])
2592  {
2593  if (mask[cells[l][1]] > 0)
2594  {
2595  if (ind == 0)
2596  j = cells[l][3];
2597  else
2598  k = cells[l][3];
2599  ind++;
2600  }
2601 
2602  if (mask[cells[l][0]] > 0)
2603  {
2604  if (ind == 0)
2605  j = cells[l][4];
2606  else
2607  k = cells[l][4];
2608  ind++;
2609  }
2610  }
2611 
2612  if (i == cells[l][3])
2613  {
2614  j = cells[l][1];
2615  k = cells[l][2];
2616  }
2617 
2618  if (i == cells[l][4])
2619  {
2620  j = cells[l][2];
2621  k = cells[l][0];
2622  }
2623 
2624  if (i == cells[l][5])
2625  {
2626  j = cells[l][0];
2627  k = cells[l][1];
2628  }
2629  break;
2630 
2631  default:
2632  libmesh_error_msg("Unrecognized nvert = " << nvert);
2633  }
2634  } // end boundary connectivity
2635 
2636  // lines
2637  Real del1 = R[j][0] - R[i][0];
2638  Real del2 = R[i][0] - R[k][0];
2639 
2640  if ((std::abs(del1) < eps) &&
2641  (std::abs(del2) < eps))
2642  {
2643  constr[I*4] = 1;
2644  constr[I*4+1] = 0;
2645  constr[I*4+2] = R[i][0];
2646  constr[I*4+3] = R[i][1];
2647  }
2648  else
2649  {
2650  del1 = R[j][1] - R[i][1];
2651  del2 = R[i][1] - R[k][1];
2652  if ((std::abs(del1) < eps) &&
2653  (std::abs(del2) < eps))
2654  {
2655  constr[I*4] = 0;
2656  constr[I*4+1] = 1;
2657  constr[I*4+2] = R[i][0];
2658  constr[I*4+3] = R[i][1];
2659  }
2660  else
2661  {
2662  J =
2663  (R[i][0]-R[j][0])*(R[k][1]-R[j][1]) -
2664  (R[k][0]-R[j][0])*(R[i][1]-R[j][1]);
2665 
2666  if (std::abs(J) < eps)
2667  {
2668  constr[I*4] = 1./(R[k][0]-R[j][0]);
2669  constr[I*4+1] = -1./(R[k][1]-R[j][1]);
2670  constr[I*4+2] = R[i][0];
2671  constr[I*4+3] = R[i][1];
2672  }
2673  else
2674  {
2675  // circle
2676  Real x0 = ((R[k][1]-R[j][1])*(R[i][0]*R[i][0]-R[j][0]*R[j][0]+R[i][1]*R[i][1]-R[j][1]*R[j][1]) +
2677  (R[j][1]-R[i][1])*(R[k][0]*R[k][0]-R[j][0]*R[j][0]+R[k][1]*R[k][1]-R[j][1]*R[j][1]))/(2*J);
2678  Real y0 = ((R[j][0]-R[k][0])*(R[i][0]*R[i][0]-R[j][0]*R[j][0]+R[i][1]*R[i][1]-R[j][1]*R[j][1]) +
2679  (R[i][0]-R[j][0])*(R[k][0]*R[k][0]-R[j][0]*R[j][0]+R[k][1]*R[k][1]-R[j][1]*R[j][1]))/(2*J);
2680  constr[I*4] = x0;
2681  constr[I*4+1] = y0;
2682  constr[I*4+2] = (R[i][0]-x0)*(R[i][0]-x0) + (R[i][1]-y0)*(R[i][1]-y0);
2683  }
2684  }
2685  }
2686  }
2687 
2688  // for (int i=0; i<NCN; i++){
2689  // for (int j=0; j<4; j++) fprintf(stdout," %e ",constr[4*i+j]);
2690  // fprintf(stdout, "constr %d node %d \n",i,Bind[i]);
2691  // }
2692 
2693  // initial values
2694  for (int i=0; i<NCN; i++)
2695  lam[i] = 0;
2696 
2697  // Eventual return value
2698  Real nonzero = 0.;
2699 
2700  // Temporary result variable
2701  Real qq = 0.;
2702 
2703  for (int nz=0; nz<5; nz++)
2704  {
2705  // find H and -grad J
2706  nonzero = 0.;
2707  Jpr = 0;
2708  for (dof_id_type i=0; i<2*_n_nodes; i++)
2709  {
2710  b[i] = 0;
2711  hm[i] = 0;
2712  }
2713 
2714  for (dof_id_type i=0; i<_n_cells; i++)
2715  {
2716  int nvert = 0;
2717 
2718  while (cells[i][nvert] >= 0)
2719  nvert++;
2720 
2721  for (int j=0; j<nvert; j++)
2722  {
2723  G[i][j] = 0;
2724  if (adp < 0)
2725  for (auto k : IntRange<int>(0, std::abs(adp)))
2726  G[i][j] += afun[i*(-adp) + k];
2727  }
2728 
2729  for (int index=0; index<2; index++)
2730  for (int k=0; k<nvert; k++)
2731  {
2732  F[index][k] = 0;
2733  for (int j=0; j<nvert; j++)
2734  W[index][k][j] = 0;
2735  }
2736 
2737  if (mcells[i] >= 0)
2738  Jpr += localP(W, F, R, cells[i], mask, epsilon, w, nvert, H[i],
2739  me, vol, 0, lVmin, lqmin, adp, afun, G[i]);
2740 
2741  else
2742  {
2743  for (unsigned index=0; index<2; index++)
2744  for (int j=0; j<nvert; j++)
2745  W[index][j][j] = 1;
2746  }
2747 
2748  for (unsigned index=0; index<2; index++)
2749  for (int l=0; l<nvert; l++)
2750  {
2751  // diagonal Hessian
2752  hm[cells[i][l] + index*_n_nodes] += W[index][l][l];
2753  b[cells[i][l] + index*_n_nodes] -= F[index][l];
2754  }
2755  }
2756 
2757  // ||grad J||_2
2758  for (dof_id_type i=0; i<2*_n_nodes; i++)
2759  nonzero += b[i]*b[i];
2760 
2761  // solve for Plam
2762  for (int I=0; I<NCN; I++)
2763  {
2764  int i = Bind[I];
2765  Real
2766  Bx = 0.,
2767  By = 0.,
2768  g = 0.;
2769 
2770  if (constr[4*I+3] < 0.5/eps)
2771  {
2772  Bx = constr[4*I];
2773  By = constr[4*I+1];
2774  g = (R[i][0]-constr[4*I+2])*constr[4*I] + (R[i][1]-constr[4*I+3])*constr[4*I+1];
2775  }
2776  else
2777  {
2778  Bx = 2*(R[i][0] - constr[4*I]);
2779  By = 2*(R[i][1] - constr[4*I+1]);
2780  hm[i] += 2*lam[I];
2781  hm[i+_n_nodes] += 2*lam[I];
2782  g =
2783  (R[i][0] - constr[4*I])*(R[i][0]-constr[4*I]) +
2784  (R[i][1] - constr[4*I+1])*(R[i][1]-constr[4*I+1]) - constr[4*I+2];
2785  }
2786 
2787  Jpr += lam[I]*g;
2788  qq = Bx*b[i]/hm[i] + By*b[i+_n_nodes]/hm[i+_n_nodes] - g;
2789  Real a = Bx*Bx/hm[i] + By*By/hm[i+_n_nodes];
2790 
2791  if (a != 0)
2792  Plam[I] = qq/a;
2793  else
2794  _logfile << "error: B^TH-1B is degenerate" << std::endl;
2795 
2796  b[i] -= Plam[I]*Bx;
2797  b[i+_n_nodes] -= Plam[I]*By;
2798  Plam[I] -= lam[I];
2799  }
2800 
2801  // solve for P
2802  for (dof_id_type i=0; i<_n_nodes; i++)
2803  {
2804  P[i][0] = b[i]/hm[i];
2805  P[i][1] = b[i+_n_nodes]/hm[i+_n_nodes];
2806  }
2807 
2808  // correct solution
2809  for (dof_id_type i=0; i<_n_nodes; i++)
2810  for (unsigned j=0; j<2; j++)
2811  if ((std::abs(P[i][j]) < eps) || (mask[i] == 1))
2812  P[i][j] = 0;
2813 
2814  // P is determined
2815  if (msglev >= 3)
2816  {
2817  for (dof_id_type i=0; i<_n_nodes; i++)
2818  for (unsigned j=0; j<2; j++)
2819  if (P[i][j] != 0)
2820  _logfile << "P[" << i << "][" << j << "]=" << P[i][j] << " ";
2821  }
2822 
2823  // local minimization problem, determination of tau
2824  if (msglev >= 3)
2825  _logfile << "dJ=" << std::sqrt(nonzero) << " L0=" << Jpr << std::endl;
2826 
2827  L = 1.e32;
2828  int j = 1;
2829 
2830  while ((Jpr <= L) && (j > -30))
2831  {
2832  j = j-1;
2833  tau = pow(2.,j);
2834  for (dof_id_type i=0; i<_n_nodes; i++)
2835  for (unsigned k=0; k<2; k++)
2836  Rpr[i][k] = R[i][k] + tau*P[i][k];
2837 
2838  J = 0;
2839  gVmin = 1.e32;
2840  gemax = -1.e32;
2841  gqmin = 1.e32;
2842  for (dof_id_type i=0; i<_n_cells; i++)
2843  if (mcells[i] >= 0)
2844  {
2845  int nvert = 0;
2846  while (cells[i][nvert] >= 0)
2847  nvert++;
2848 
2849  lemax = localP(W, F, Rpr, cells[i], mask, epsilon, w, nvert, H[i], me, vol, 1, lVmin,
2850  lqmin, adp, afun, G[i]);
2851  J += lemax;
2852 
2853  if (gVmin > lVmin)
2854  gVmin = lVmin;
2855 
2856  if (gemax < lemax)
2857  gemax = lemax;
2858 
2859  if (gqmin > lqmin)
2860  gqmin = lqmin;
2861  }
2862 
2863  L = J;
2864 
2865  // constraints contribution
2866  for (int I=0; I<NCN; I++)
2867  {
2868  int i = Bind[I];
2869  Real g = 0.;
2870 
2871  if (constr[4*I+3] < 0.5/eps)
2872  g = (Rpr[i][0] - constr[4*I+2])*constr[4*I] + (Rpr[i][1]-constr[4*I+3])*constr[4*I+1];
2873 
2874  else
2875  g = (Rpr[i][0]-constr[4*I])*(Rpr[i][0]-constr[4*I]) +
2876  (Rpr[i][1]-constr[4*I+1])*(Rpr[i][1]-constr[4*I+1]) - constr[4*I+2];
2877 
2878  L += (lam[I] + tau*Plam[I])*g;
2879  }
2880 
2881  // end of constraints
2882  if (msglev >= 3)
2883  _logfile << " tau=" << tau << " J=" << J << std::endl;
2884  } // end while
2885 
2886  if (j == -30)
2887  T = 0;
2888  else
2889  {
2890  Jpr = L;
2891  qq = J;
2892  for (dof_id_type i=0; i<_n_nodes; i++)
2893  for (unsigned k=0; k<2; k++)
2894  Rpr[i][k] = R[i][k] + tau*0.5*P[i][k];
2895 
2896  J = 0;
2897  gVmin0 = 1.e32;
2898  gemax0 = -1.e32;
2899  gqmin0 = 1.e32;
2900 
2901  for (dof_id_type i=0; i<_n_cells; i++)
2902  if (mcells[i] >= 0)
2903  {
2904  int nvert = 0;
2905  while (cells[i][nvert] >= 0)
2906  nvert++;
2907 
2908  lemax = localP(W, F, Rpr, cells[i], mask, epsilon, w, nvert, H[i], me, vol, 1, lVmin,
2909  lqmin, adp, afun, G[i]);
2910  J += lemax;
2911 
2912  if (gVmin0 > lVmin)
2913  gVmin0 = lVmin;
2914 
2915  if (gemax0 < lemax)
2916  gemax0 = lemax;
2917 
2918  if (gqmin0 > lqmin)
2919  gqmin0 = lqmin;
2920  }
2921 
2922  L = J;
2923 
2924  // constraints contribution
2925  for (int I=0; I<NCN; I++)
2926  {
2927  int i = Bind[I];
2928  Real g = 0.;
2929 
2930  if (constr[4*I+3] < 0.5/eps)
2931  g = (Rpr[i][0]-constr[4*I+2])*constr[4*I] + (Rpr[i][1]-constr[4*I+3])*constr[4*I+1];
2932 
2933  else
2934  g = (Rpr[i][0]-constr[4*I])*(Rpr[i][0]-constr[4*I]) +
2935  (Rpr[i][1]-constr[4*I+1])*(Rpr[i][1]-constr[4*I+1]) - constr[4*I+2];
2936 
2937  L += (lam[I] + tau*0.5*Plam[I])*g;
2938  }
2939  // end of constraints
2940  }
2941 
2942  if (Jpr > L)
2943  {
2944  T = 0.5*tau;
2945  // Set return values by reference
2946  Vmin = gVmin0;
2947  emax = gemax0;
2948  qmin = gqmin0;
2949  }
2950  else
2951  {
2952  T = tau;
2953  L = Jpr;
2954  J = qq;
2955  // Set return values by reference
2956  Vmin = gVmin;
2957  emax = gemax;
2958  qmin = gqmin;
2959  }
2960 
2961  for (dof_id_type i=0; i<_n_nodes; i++)
2962  for (unsigned k=0; k<2; k++)
2963  R[i][k] += T*P[i][k];
2964 
2965  for (int i=0; i<NCN; i++)
2966  lam[i] += T*Plam[i];
2967 
2968  } // end Lagrangian iter
2969 
2970  if (msglev >= 2)
2971  _logfile << "tau=" << T << ", J=" << J << ", L=" << L << std::endl;
2972 
2973  return std::sqrt(nonzero);
2974 }

References _logfile, _n_cells, _n_nodes, std::abs(), localP(), std::pow(), libMesh::Real, and std::sqrt().

Referenced by full_smooth().

◆ minq()

Real libMesh::VariationalMeshSmoother::minq ( const Array2D< Real > &  R,
const Array2D< int > &  cells,
const std::vector< int > &  mcells,
int  me,
const Array3D< Real > &  H,
Real &  vol,
Real &  Vmin 
)
private

Definition at line 1484 of file mesh_smoother_vsmoother.C.

1491 {
1492  std::vector<Real> K(9);
1493  Array2D<Real> Q(3, 3*_dim + _dim%2);
1494 
1495  Real v = 0;
1496  Real vmin = 1.e32;
1497  Real gqmin = 1.e32;
1498 
1499  for (dof_id_type ii=0; ii<_n_cells; ii++)
1500  if (mcells[ii] >= 0)
1501  {
1502  if (_dim == 2)
1503  {
1504  // 2D
1505  if (cells[ii][3] == -1)
1506  {
1507  // tri
1508  basisA(Q, 3, K, H[ii], me);
1509 
1510  std::vector<Real> a1(3), a2(3);
1511  for (int k=0; k<2; k++)
1512  for (int l=0; l<3; l++)
1513  {
1514  a1[k] += Q[k][l]*R[cells[ii][l]][0];
1515  a2[k] += Q[k][l]*R[cells[ii][l]][1];
1516  }
1517 
1518  Real det = jac2(a1[0],a1[1],a2[0],a2[1]);
1519  if (gqmin > det)
1520  gqmin = det;
1521 
1522  if (vmin > det)
1523  vmin = det;
1524 
1525  v += det;
1526  }
1527  else if (cells[ii][4] == -1)
1528  {
1529  // quad
1530  Real vcell = 0.;
1531  for (int i=0; i<2; i++)
1532  {
1533  K[0] = i;
1534  for (int j=0; j<2; j++)
1535  {
1536  K[1] = j;
1537  basisA(Q, 4, K, H[ii], me);
1538 
1539  std::vector<Real> a1(3), a2(3);
1540  for (int k=0; k<2; k++)
1541  for (int l=0; l<4; l++)
1542  {
1543  a1[k] += Q[k][l]*R[cells[ii][l]][0];
1544  a2[k] += Q[k][l]*R[cells[ii][l]][1];
1545  }
1546 
1547  Real det = jac2(a1[0],a1[1],a2[0],a2[1]);
1548  if (gqmin > det)
1549  gqmin = det;
1550 
1551  v += 0.25*det;
1552  vcell += 0.25*det;
1553  }
1554  }
1555  if (vmin > vcell)
1556  vmin = vcell;
1557  }
1558  else
1559  {
1560  // quad tri
1561  Real vcell = 0.;
1562  for (int i=0; i<3; i++)
1563  {
1564  K[0] = i*0.5;
1565  int k = i/2;
1566  K[1] = static_cast<Real>(k);
1567 
1568  for (int j=0; j<3; j++)
1569  {
1570  K[2] = j*0.5;
1571  k = j/2;
1572  K[3] = static_cast<Real>(k);
1573  basisA(Q, 6, K, H[ii], me);
1574 
1575  std::vector<Real> a1(3), a2(3);
1576  for (int k2=0; k2<2; k2++)
1577  for (int l=0; l<6; l++)
1578  {
1579  a1[k2] += Q[k2][l]*R[cells[ii][l]][0];
1580  a2[k2] += Q[k2][l]*R[cells[ii][l]][1];
1581  }
1582 
1583  Real det = jac2(a1[0], a1[1], a2[0], a2[1]);
1584  if (gqmin > det)
1585  gqmin = det;
1586 
1587  Real sigma = 1./24.;
1588  if (i == j)
1589  sigma = 1./12.;
1590 
1591  v += sigma*det;
1592  vcell += sigma*det;
1593  }
1594  }
1595  if (vmin > vcell)
1596  vmin = vcell;
1597  }
1598  }
1599  if (_dim == 3)
1600  {
1601  // 3D
1602  if (cells[ii][4] == -1)
1603  {
1604  // tetr
1605  basisA(Q, 4, K, H[ii], me);
1606 
1607  std::vector<Real> a1(3), a2(3), a3(3);
1608  for (int k=0; k<3; k++)
1609  for (int l=0; l<4; l++)
1610  {
1611  a1[k] += Q[k][l]*R[cells[ii][l]][0];
1612  a2[k] += Q[k][l]*R[cells[ii][l]][1];
1613  a3[k] += Q[k][l]*R[cells[ii][l]][2];
1614  }
1615 
1616  Real det = jac3(a1[0], a1[1], a1[2],
1617  a2[0], a2[1], a2[2],
1618  a3[0], a3[1], a3[2]);
1619 
1620  if (gqmin > det)
1621  gqmin = det;
1622 
1623  if (vmin > det)
1624  vmin = det;
1625  v += det;
1626  }
1627  else if (cells[ii][6] == -1)
1628  {
1629  // prism
1630  Real vcell = 0.;
1631  for (int i=0; i<2; i++)
1632  {
1633  K[0] = i;
1634  for (int j=0; j<2; j++)
1635  {
1636  K[1] = j;
1637 
1638  for (int k=0; k<3; k++)
1639  {
1640  K[2] = 0.5*k;
1641  K[3] = static_cast<Real>(k%2);
1642  basisA(Q, 6, K, H[ii], me);
1643 
1644  std::vector<Real> a1(3), a2(3), a3(3);
1645  for (int kk=0; kk<3; kk++)
1646  for (int ll=0; ll<6; ll++)
1647  {
1648  a1[kk] += Q[kk][ll]*R[cells[ii][ll]][0];
1649  a2[kk] += Q[kk][ll]*R[cells[ii][ll]][1];
1650  a3[kk] += Q[kk][ll]*R[cells[ii][ll]][2];
1651  }
1652 
1653  Real det = jac3(a1[0], a1[1], a1[2],
1654  a2[0], a2[1], a2[2],
1655  a3[0], a3[1], a3[2]);
1656  if (gqmin > det)
1657  gqmin = det;
1658 
1659  Real sigma = 1./12.;
1660  v += sigma*det;
1661  vcell += sigma*det;
1662  }
1663  }
1664  }
1665  if (vmin > vcell)
1666  vmin = vcell;
1667  }
1668  else if (cells[ii][8] == -1)
1669  {
1670  // hex
1671  Real vcell = 0.;
1672  for (int i=0; i<2; i++)
1673  {
1674  K[0] = i;
1675  for (int j=0; j<2; j++)
1676  {
1677  K[1] = j;
1678  for (int k=0; k<2; k++)
1679  {
1680  K[2] = k;
1681  for (int l=0; l<2; l++)
1682  {
1683  K[3] = l;
1684  for (int m=0; m<2; m++)
1685  {
1686  K[4] = m;
1687  for (int nn=0; nn<2; nn++)
1688  {
1689  K[5] = nn;
1690  basisA(Q, 8, K, H[ii], me);
1691 
1692  std::vector<Real> a1(3), a2(3), a3(3);
1693  for (int kk=0; kk<3; kk++)
1694  for (int ll=0; ll<8; ll++)
1695  {
1696  a1[kk] += Q[kk][ll]*R[cells[ii][ll]][0];
1697  a2[kk] += Q[kk][ll]*R[cells[ii][ll]][1];
1698  a3[kk] += Q[kk][ll]*R[cells[ii][ll]][2];
1699  }
1700 
1701  Real det = jac3(a1[0], a1[1], a1[2],
1702  a2[0], a2[1], a2[2],
1703  a3[0], a3[1], a3[2]);
1704 
1705  if (gqmin > det)
1706  gqmin = det;
1707 
1708  Real sigma = 0.;
1709 
1710  if ((i==nn) && (j==l) && (k==m))
1711  sigma = 1./27.;
1712 
1713  if (((i==nn) && (j==l) && (k!=m)) ||
1714  ((i==nn) && (j!=l) && (k==m)) ||
1715  ((i!=nn) && (j==l) && (k==m)))
1716  sigma = 1./54.;
1717 
1718  if (((i==nn) && (j!=l) && (k!=m)) ||
1719  ((i!=nn) && (j!=l) && (k==m)) ||
1720  ((i!=nn) && (j==l) && (k!=m)))
1721  sigma = 1./108.;
1722 
1723  if ((i!=nn) && (j!=l) && (k!=m))
1724  sigma = 1./216.;
1725 
1726  v += sigma*det;
1727  vcell += sigma*det;
1728  }
1729  }
1730  }
1731  }
1732  }
1733  }
1734 
1735  if (vmin > vcell)
1736  vmin = vcell;
1737  }
1738  else
1739  {
1740  // quad tetr
1741  Real vcell = 0.;
1742  for (int i=0; i<4; i++)
1743  {
1744  for (int j=0; j<4; j++)
1745  {
1746  for (int k=0; k<4; k++)
1747  {
1748  switch (i)
1749  {
1750  case 0:
1751  K[0] = 0;
1752  K[1] = 0;
1753  K[2] = 0;
1754  break;
1755 
1756  case 1:
1757  K[0] = 1;
1758  K[1] = 0;
1759  K[2] = 0;
1760  break;
1761 
1762  case 2:
1763  K[0] = 0.5;
1764  K[1] = 1;
1765  K[2] = 0;
1766  break;
1767 
1768  case 3:
1769  K[0] = 0.5;
1770  K[1] = 1./3.;
1771  K[2] = 1;
1772  break;
1773 
1774  default:
1775  break;
1776  }
1777  switch (j)
1778  {
1779  case 0:
1780  K[3] = 0;
1781  K[4] = 0;
1782  K[5] = 0;
1783  break;
1784 
1785  case 1:
1786  K[3] = 1;
1787  K[4] = 0;
1788  K[5] = 0;
1789  break;
1790 
1791  case 2:
1792  K[3] = 0.5;
1793  K[4] = 1;
1794  K[5] = 0;
1795  break;
1796 
1797  case 3:
1798  K[3] = 0.5;
1799  K[4] = 1./3.;
1800  K[5] = 1;
1801  break;
1802 
1803  default:
1804  break;
1805  }
1806  switch (k)
1807  {
1808  case 0:
1809  K[6] = 0;
1810  K[7] = 0;
1811  K[8] = 0;
1812  break;
1813 
1814  case 1:
1815  K[6] = 1;
1816  K[7] = 0;
1817  K[8] = 0;
1818  break;
1819 
1820  case 2:
1821  K[6] = 0.5;
1822  K[7] = 1;
1823  K[8] = 0;
1824  break;
1825 
1826  case 3:
1827  K[6] = 0.5;
1828  K[7] = 1./3.;
1829  K[8] = 1;
1830  break;
1831 
1832  default:
1833  break;
1834  }
1835 
1836  basisA(Q, 10, K, H[ii], me);
1837 
1838  std::vector<Real> a1(3), a2(3), a3(3);
1839  for (int kk=0; kk<3; kk++)
1840  for (int ll=0; ll<10; ll++)
1841  {
1842  a1[kk] += Q[kk][ll]*R[cells[ii][ll]][0];
1843  a2[kk] += Q[kk][ll]*R[cells[ii][ll]][1];
1844  a3[kk] += Q[kk][ll]*R[cells[ii][ll]][2];
1845  }
1846 
1847  Real det = jac3(a1[0], a1[1], a1[2],
1848  a2[0], a2[1], a2[2],
1849  a3[0], a3[1], a3[2]);
1850  if (gqmin > det)
1851  gqmin = det;
1852 
1853  Real sigma = 0.;
1854 
1855  if ((i==j) && (j==k))
1856  sigma = 1./120.;
1857 
1858  else if ((i==j) || (j==k) || (i==k))
1859  sigma = 1./360.;
1860 
1861  else
1862  sigma = 1./720.;
1863 
1864  v += sigma*det;
1865  vcell += sigma*det;
1866  }
1867  }
1868  }
1869  if (vmin > vcell)
1870  vmin = vcell;
1871  }
1872  }
1873  }
1874 
1875  // Fill in return value references
1876  vol = v/static_cast<Real>(_n_cells);
1877  Vmin = vmin;
1878 
1879  return gqmin;
1880 }

References _dim, _n_cells, basisA(), jac2(), jac3(), and libMesh::Real.

Referenced by full_smooth().

◆ pcg_ic0()

int libMesh::VariationalMeshSmoother::pcg_ic0 ( int  n,
const std::vector< int > &  ia,
const std::vector< int > &  ja,
const std::vector< Real > &  a,
const std::vector< Real > &  u,
std::vector< Real > &  x,
const std::vector< Real > &  b,
std::vector< Real > &  r,
std::vector< Real > &  p,
std::vector< Real > &  z,
Real  eps,
int  maxite,
int  msglev 
)
private

Definition at line 3806 of file mesh_smoother_vsmoother.C.

3819 {
3820  // Return value
3821  int i = 0;
3822 
3823  Real
3824  rhr = 0.,
3825  rhr0 = 0.,
3826  rhrold = 0.,
3827  rr0 = 0.;
3828 
3829  for (i=0; i<=maxite; i++)
3830  {
3831  if (i == 0)
3832  p = x;
3833 
3834  // z=Ap
3835  for (int ii=0; ii<n; ii++)
3836  z[ii] = 0.;
3837 
3838  for (int ii=0; ii<n; ii++)
3839  {
3840  z[ii] += a[ia[ii]]*p[ii];
3841 
3842  for (int j=ia[ii]+1; j<ia[ii+1]; j++)
3843  {
3844  z[ii] += a[j]*p[ja[j]-1];
3845  z[ja[j]-1] += a[j]*p[ii];
3846  }
3847  }
3848 
3849  if (i == 0)
3850  for (int k=0; k<n; k++)
3851  r[k] = b[k] - z[k]; // r(0) = b - Ax(0)
3852 
3853  if (i > 0)
3854  {
3855  Real pap = 0.;
3856  for (int k=0; k<n; k++)
3857  pap += p[k]*z[k];
3858 
3859  Real alpha = rhr/pap;
3860  for (int k=0; k<n; k++)
3861  {
3862  x[k] += p[k]*alpha;
3863  r[k] -= z[k]*alpha;
3864  }
3865  rhrold = rhr;
3866  }
3867 
3868  // z = H * r
3869  for (int ii=0; ii<n; ii++)
3870  z[ii] = r[ii]*u[ii];
3871 
3872  if (i == 0)
3873  p = z;
3874 
3875  rhr = 0.;
3876  Real rr = 0.;
3877  for (int k=0; k<n; k++)
3878  {
3879  rhr += r[k]*z[k];
3880  rr += r[k]*r[k];
3881  }
3882 
3883  if (i == 0)
3884  {
3885  rhr0 = rhr;
3886  rr0 = rr;
3887  }
3888 
3889  if (msglev > 0)
3890  _logfile << i
3891  << " ) rHr ="
3892  << std::sqrt(rhr/rhr0)
3893  << " rr ="
3894  << std::sqrt(rr/rr0)
3895  << std::endl;
3896 
3897  if (rr <= eps*eps*rr0)
3898  return i;
3899 
3900  if (i >= maxite)
3901  return i;
3902 
3903  if (i > 0)
3904  {
3905  Real beta = rhr/rhrold;
3906  for (int k=0; k<n; k++)
3907  p[k] = z[k] + p[k]*beta;
3908  }
3909  } // end for i
3910 
3911  return i;
3912 }

References _logfile, libMesh::Real, and std::sqrt().

Referenced by solver().

◆ pcg_par_check()

int libMesh::VariationalMeshSmoother::pcg_par_check ( int  n,
const std::vector< int > &  ia,
const std::vector< int > &  ja,
const std::vector< Real > &  a,
Real  eps,
int  maxite,
int  msglev 
)
private

Definition at line 3664 of file mesh_smoother_vsmoother.C.

3671 {
3672  int i, j, jj, k;
3673 
3674  if (n <= 0)
3675  {
3676  if (msglev > 0)
3677  _logfile << "sol_pcg: incorrect input parameter: n =" << n << std::endl;
3678  return -1;
3679  }
3680 
3681  if (ia[0] != 0)
3682  {
3683  if (msglev > 0)
3684  _logfile << "sol_pcg: incorrect input parameter: ia(1) =" << ia[0] << std::endl;
3685  return -2;
3686  }
3687 
3688  for (i=0; i<n; i++)
3689  {
3690  if (ia[i+1] < ia[i])
3691  {
3692  if (msglev >= 1)
3693  _logfile << "sol_pcg: incorrect input parameter: i ="
3694  << i+1
3695  << "; ia(i) ="
3696  << ia[i]
3697  << " must be <= a(i+1) ="
3698  << ia[i+1]
3699  << std::endl;
3700  return -2;
3701  }
3702  }
3703 
3704  for (i=0; i<n; i++)
3705  {
3706  if (ja[ia[i]] != (i+1))
3707  {
3708  if (msglev >= 1)
3709  _logfile << "sol_pcg: incorrect input parameter: i ="
3710  << i+1
3711  << " ; ia(i) ="
3712  << ia[i]
3713  << " ; absence of diagonal entry; k ="
3714  << ia[i]+1
3715  << " ; the value ja(k) ="
3716  << ja[ia[i]]
3717  << " must be equal to i"
3718  << std::endl;
3719 
3720  return -3;
3721  }
3722 
3723  jj = 0;
3724  for (k=ia[i]; k<ia[i+1]; k++)
3725  {
3726  j=ja[k];
3727  if ((j<(i+1)) || (j>n))
3728  {
3729  if (msglev >= 1)
3730  _logfile << "sol_pcg: incorrect input parameter: i ="
3731  << i+1
3732  << " ; ia(i) ="
3733  << ia[i]
3734  << " ; a(i+1) ="
3735  << ia[i+1]
3736  << " ; k ="
3737  << k+1
3738  << " ; the value ja(k) ="
3739  << ja[k]
3740  << " is out of range"
3741  << std::endl;
3742 
3743  return -3;
3744  }
3745  if (j <= jj)
3746  {
3747  if (msglev >= 1)
3748  _logfile << "sol_pcg: incorrect input parameter: i ="
3749  << i+1
3750  << " ; ia(i) ="
3751  << ia[i]
3752  << " ; a(i+1) ="
3753  << ia[i+1]
3754  << " ; k ="
3755  << k+1
3756  << " ; the value ja(k) ="
3757  << ja[k]
3758  << " must be less than ja(k+1) ="
3759  << ja[k+1]
3760  << std::endl;
3761 
3762  return -3;
3763  }
3764  jj = j;
3765  }
3766  }
3767 
3768  for (i=0; i<n; i++)
3769  {
3770  if (a[ia[i]] <= 0.)
3771  {
3772  if (msglev >= 1)
3773  _logfile << "sol_pcg: incorrect input parameter: i ="
3774  << i+1
3775  << " ; ia(i) ="
3776  << ia[i]
3777  << " ; the diagonal entry a(ia(i)) ="
3778  << a[ia[i]]
3779  << " must be > 0"
3780  << std::endl;
3781 
3782  return -4;
3783  }
3784  }
3785 
3786  if (eps < 0.)
3787  {
3788  if (msglev > 0)
3789  _logfile << "sol_pcg: incorrect input parameter: eps =" << eps << std::endl;
3790  return -7;
3791  }
3792 
3793  if (maxite < 1)
3794  {
3795  if (msglev > 0)
3796  _logfile << "sol_pcg: incorrect input parameter: maxite =" << maxite << std::endl;
3797  return -8;
3798  }
3799 
3800  return 0;
3801 }

References _logfile.

Referenced by solver().

◆ read_adp()

int libMesh::VariationalMeshSmoother::read_adp ( std::vector< Real > &  afun)
private

Definition at line 573 of file mesh_smoother_vsmoother.C.

574 {
575  std::vector<float> & adapt_data = *_adapt_data;
576 
577  if (_area_of_interest)
578  adjust_adapt_data();
579 
580  std::size_t m = adapt_data.size();
581 
582  std::size_t j =0 ;
583  for (std::size_t i=0; i<m; i++)
584  if (adapt_data[i] != 0)
585  {
586  afun[j] = adapt_data[i];
587  j++;
588  }
589 
590  return 0;
591 }

References _adapt_data, _area_of_interest, and adjust_adapt_data().

Referenced by full_smooth().

◆ readgr()

int libMesh::VariationalMeshSmoother::readgr ( Array2D< Real > &  R,
std::vector< int > &  mask,
Array2D< int > &  cells,
std::vector< int > &  mcells,
std::vector< int > &  edges,
std::vector< int > &  hnodes 
)
private

Definition at line 269 of file mesh_smoother_vsmoother.C.

275 {
276  libMesh::out << "Starting readgr" << std::endl;
277  // add error messages where format can be inconsistent
278 
279  // Create a quickly-searchable list of boundary nodes.
280  std::unordered_set<dof_id_type> boundary_node_ids =
281  MeshTools::find_boundary_nodes (_mesh);
282 
283  // Grab node coordinates and set mask
284  {
285  // Only compute the node to elem map once
286  std::unordered_map<dof_id_type, std::vector<const Elem *>> nodes_to_elem_map;
287  MeshTools::build_nodes_to_elem_map(_mesh, nodes_to_elem_map);
288 
289  int i = 0;
290  for (auto & node : _mesh.node_ptr_range())
291  {
292  // Get a reference to the node
293  Node & node_ref = *node;
294 
295  // For each node grab its X Y [Z] coordinates
296  for (unsigned int j=0; j<_dim; j++)
297  R[i][j] = node_ref(j);
298 
299  // Set the Proper Mask
300  // Internal nodes are 0
301  // Immovable boundary nodes are 1
302  // Movable boundary nodes are 2
303  if (boundary_node_ids.count(i))
304  {
305  // Only look for sliding edge nodes in 2D
306  if (_dim == 2)
307  {
308  // Find all the nodal neighbors... that is the nodes directly connected
309  // to this node through one edge
310  std::vector<const Node *> neighbors;
311  MeshTools::find_nodal_neighbors(_mesh, node_ref, nodes_to_elem_map, neighbors);
312 
313  // Grab the x,y coordinates
314  Real x = node_ref(0);
315  Real y = node_ref(1);
316 
317  // Theta will represent the atan2 angle (meaning with the proper quadrant in mind)
318  // of the neighbor node in a system where the current node is at the origin
319  Real theta = 0;
320  std::vector<Real> thetas;
321 
322  // Calculate the thetas
323  for (const auto & neighbor : neighbors)
324  {
325  // Note that the x and y values of this node are subtracted off
326  // this centers the system around this node
327  theta = atan2((*neighbor)(1)-y, (*neighbor)(0)-x);
328 
329  // Save it for later
330  thetas.push_back(theta);
331  }
332 
333  // Assume the node is immovable... then prove otherwise
334  mask[i] = 1;
335 
336  // Search through neighbor nodes looking for two that form a straight line with this node
337  for (auto a : IntRange<std::size_t>(0, thetas.size()-1))
338  {
339  // Only try each pairing once
340  for (auto b : IntRange<std::size_t>(a+1, thetas.size()))
341  {
342  // Find if the two neighbor nodes angles are 180 degrees (pi) off of each other (withing a tolerance)
343  // In order to make this a true movable boundary node... the two that forma straight line with
344  // it must also be on the boundary
345  if (boundary_node_ids.count(neighbors[a]->id()) &&
346  boundary_node_ids.count(neighbors[b]->id()) &&
347  (
348  (std::abs(thetas[a] - (thetas[b] + (libMesh::pi))) < .001) ||
349  (std::abs(thetas[a] - (thetas[b] - (libMesh::pi))) < .001)
350  )
351  )
352  {
353  // if ((*(*it))(1) > 0.25 || (*(*it))(0) > .7 || (*(*it))(0) < .01)
354  mask[i] = 2;
355  }
356 
357  }
358  }
359  }
360  else // In 3D set all boundary nodes to be fixed
361  mask[i] = 1;
362  }
363  else
364  mask[i] = 0; // Internal Node
365 
366  // libMesh::out << "Node: " << i << " Mask: " << mask[i] << std::endl;
367  i++;
368  }
369  }
370 
371  // Grab the connectivity
372  // FIXME: Generalize this!
373  {
374  int i = 0;
375  for (const auto & elem : _mesh.active_element_ptr_range())
376  {
377  // Keep track of the number of nodes
378  // there must be 6 for 2D
379  // 10 for 3D
380  // If there are actually less than that -1 must be used
381  // to fill out the rest
382  int num = 0;
383  /*
384  int num_necessary = 0;
385 
386  if (_dim == 2)
387  num_necessary = 6;
388  else
389  num_necessary = 10;
390  */
391 
392  if (_dim == 2)
393  {
394  switch (elem->n_vertices())
395  {
396  // Grab nodes that do exist
397  case 3: // Tri
398  for (auto k : IntRange<unsigned int>(0, elem->n_vertices()))
399  cells[i][k] = elem->node_id(k);
400 
401  num = elem->n_vertices();
402  break;
403 
404  case 4: // Quad 4
405  cells[i][0] = elem->node_id(0);
406  cells[i][1] = elem->node_id(1);
407  cells[i][2] = elem->node_id(3); // Note that 2 and 3 are switched!
408  cells[i][3] = elem->node_id(2);
409  num = 4;
410  break;
411 
412  default:
413  libmesh_error_msg("Unknown number of vertices = " << elem->n_vertices());
414  }
415  }
416  else
417  {
418  // Grab nodes that do exist
419  switch (elem->n_vertices())
420  {
421  // Tet 4
422  case 4:
423  for (auto k : IntRange<unsigned int>(0, elem->n_vertices()))
424  cells[i][k] = elem->node_id(k);
425  num = elem->n_vertices();
426  break;
427 
428  // Hex 8
429  case 8:
430  cells[i][0] = elem->node_id(0);
431  cells[i][1] = elem->node_id(1);
432  cells[i][2] = elem->node_id(3); // Note that 2 and 3 are switched!
433  cells[i][3] = elem->node_id(2);
434 
435  cells[i][4] = elem->node_id(4);
436  cells[i][5] = elem->node_id(5);
437  cells[i][6] = elem->node_id(7); // Note that 6 and 7 are switched!
438  cells[i][7] = elem->node_id(6);
439  num=8;
440  break;
441 
442  default:
443  libmesh_error_msg("Unknown 3D element with = " << elem->n_vertices() << " vertices.");
444  }
445  }
446 
447  // Fill in the rest with -1
448  for (auto j : IntRange<int>(num, cast_int<int>(cells[i].size())))
449  cells[i][j] = -1;
450 
451  // Mask it with 0 to state that this is an active element
452  // FIXME: Could be something other than zero
453  mcells[i] = 0;
454  i++;
455  }
456  }
457 
458  // Grab hanging node connectivity
459  {
460  std::map<dof_id_type, std::vector<dof_id_type>>::iterator
461  it = _hanging_nodes.begin(),
462  end = _hanging_nodes.end();
463 
464  for (int i=0; it!=end; it++)
465  {
466  libMesh::out << "Parent 1: " << (it->second)[0] << std::endl;
467  libMesh::out << "Parent 2: " << (it->second)[1] << std::endl;
468  libMesh::out << "Hanging Node: " << it->first << std::endl << std::endl;
469 
470  // First Parent
471  edges[2*i] = (it->second)[1];
472 
473  // Second Parent
474  edges[2*i+1] = (it->second)[0];
475 
476  // Hanging Node
477  hnodes[i] = it->first;
478 
479  i++;
480  }
481  }
482  libMesh::out << "Finished readgr" << std::endl;
483 
484  return 0;
485 }

References _dim, _hanging_nodes, libMesh::MeshSmoother::_mesh, std::abs(), libMesh::MeshBase::active_element_ptr_range(), libMesh::MeshTools::build_nodes_to_elem_map(), end, libMesh::MeshTools::find_boundary_nodes(), libMesh::MeshTools::find_nodal_neighbors(), libMesh::MeshBase::node_ptr_range(), libMesh::out, libMesh::pi, and libMesh::Real.

Referenced by metr_data_gen(), and smooth().

◆ readmetr()

int libMesh::VariationalMeshSmoother::readmetr ( std::string  name,
Array3D< Real > &  H 
)
private

Definition at line 490 of file mesh_smoother_vsmoother.C.

492 {
493  std::ifstream infile(name.c_str());
494  std::string dummy;
495 
496  for (dof_id_type i=0; i<_n_cells; i++)
497  for (unsigned j=0; j<_dim; j++)
498  {
499  for (unsigned k=0; k<_dim; k++)
500  infile >> H[i][j][k];
501 
502  // Read to end of line and discard
503  std::getline(infile, dummy);
504  }
505 
506  return 0;
507 }

References _dim, _n_cells, and libMesh::Quality::name().

Referenced by smooth().

◆ set_generate_data()

void libMesh::VariationalMeshSmoother::set_generate_data ( bool  b)
inline

Allow user to control whether the metric is generated from the initial mesh.

Definition at line 142 of file mesh_smoother_vsmoother.h.

142 { _generate_data = b; }

References _generate_data.

◆ set_metric()

void libMesh::VariationalMeshSmoother::set_metric ( MetricType  t)
inline

Allow user to control the smoothing metric used.

Definition at line 147 of file mesh_smoother_vsmoother.h.

147 { _metric = t; }

References _metric.

◆ smooth() [1/2]

virtual void libMesh::VariationalMeshSmoother::smooth ( )
inlineoverridevirtual

Redefinition of the smooth function from the base class.

All this does is call the smooth function in this class which takes an int, using a default value of 1.

Implements libMesh::MeshSmoother.

Definition at line 125 of file mesh_smoother_vsmoother.h.

125 { _distance = this->smooth(1); }

References _distance, and smooth().

Referenced by smooth().

◆ smooth() [2/2]

Real libMesh::VariationalMeshSmoother::smooth ( unsigned int  n_iterations)

The actual smoothing function, gets called whenever the user specifies an actual number of smoothing iterations.

Definition at line 126 of file mesh_smoother_vsmoother.C.

127 {
128  // If the log file is already open, for example on subsequent calls
129  // to smooth() on the same object, we'll just keep writing to it,
130  // otherwise we'll open it...
131  if (!_logfile.is_open())
132  _logfile.open("smoother.out");
133 
134  int
135  me = _metric,
136  gr = _generate_data ? 0 : 1,
137  adp = _adaptive_func,
138  miniter = _miniter,
139  maxiter = _maxiter,
140  miniterBC = _miniterBC;
141 
142  Real theta = _theta;
143 
144  // Metric file name
145  std::string metric_filename = "smoother.metric";
146  if (gr == 0 && me > 1)
147  {
148  // grid filename
149  std::string grid_filename = "smoother.grid";
150 
151  // generate metric from initial mesh (me = 2,3)
152  metr_data_gen(grid_filename, metric_filename, me);
153  }
154 
155  // Initialize the _n_nodes and _n_cells member variables
156  this->_n_nodes = _mesh.n_nodes();
157  this->_n_cells = _mesh.n_active_elem();
158 
159  // Initialize the _n_hanging_edges member variable
160  MeshTools::find_hanging_nodes_and_parents(_mesh, _hanging_nodes);
161  this->_n_hanging_edges =
162  cast_int<dof_id_type>(_hanging_nodes.size());
163 
164  std::vector<int>
165  mask(_n_nodes),
166  edges(2*_n_hanging_edges),
167  mcells(_n_cells),
168  hnodes(_n_hanging_edges);
169 
170  Array2D<Real> R(_n_nodes, _dim);
171  Array2D<int> cells(_n_cells, 3*_dim + _dim%2);
172  Array3D<Real> H(_n_cells, _dim, _dim);
173 
174  // initial grid
175  int vms_err = readgr(R, mask, cells, mcells, edges, hnodes);
176  if (vms_err < 0)
177  {
178  _logfile << "Error reading input mesh file" << std::endl;
179  return _dist_norm;
180  }
181 
182  if (me > 1)
183  vms_err = readmetr(metric_filename, H);
184 
185  if (vms_err < 0)
186  {
187  _logfile << "Error reading metric file" << std::endl;
188  return _dist_norm;
189  }
190 
191  std::vector<int> iter(4);
192  iter[0] = miniter;
193  iter[1] = maxiter;
194  iter[2] = miniterBC;
195 
196  // grid optimization
197  _logfile << "Starting Grid Optimization" << std::endl;
198  clock_t ticks1 = clock();
199  full_smooth(R, mask, cells, mcells, edges, hnodes, theta, iter, me, H, adp, gr);
200  clock_t ticks2 = clock();
201  _logfile << "full_smooth took ("
202  << ticks2
203  << "-"
204  << ticks1
205  << ")/"
206  << CLOCKS_PER_SEC
207  << " = "
208  << static_cast<Real>(ticks2-ticks1)/static_cast<Real>(CLOCKS_PER_SEC)
209  << " seconds"
210  << std::endl;
211 
212  // save result
213  _logfile << "Saving Result" << std::endl;
214  writegr(R);
215 
216  libmesh_assert_greater (_dist_norm, 0);
217  return _dist_norm;
218 }

References _adaptive_func, _dim, _dist_norm, _generate_data, _hanging_nodes, _logfile, _maxiter, libMesh::MeshSmoother::_mesh, _metric, _miniter, _miniterBC, _n_cells, _n_hanging_edges, _n_nodes, _theta, libMesh::MeshTools::find_hanging_nodes_and_parents(), full_smooth(), metr_data_gen(), libMesh::MeshBase::n_active_elem(), libMesh::MeshBase::n_nodes(), readgr(), readmetr(), libMesh::Real, and writegr().

◆ solver()

int libMesh::VariationalMeshSmoother::solver ( int  n,
const std::vector< int > &  ia,
const std::vector< int > &  ja,
const std::vector< Real > &  a,
std::vector< Real > &  x,
const std::vector< Real > &  b,
Real  eps,
int  maxite,
int  msglev 
)
private

Definition at line 3627 of file mesh_smoother_vsmoother.C.

3636 {
3637  _logfile << "Beginning Solve " << n << std::endl;
3638 
3639  // Check parameters
3640  int info = pcg_par_check(n, ia, ja, a, eps, maxite, msglev);
3641  if (info != 0)
3642  return info;
3643 
3644  // PJ preconditioner construction
3645  std::vector<Real> u(n);
3646  for (int i=0; i<n; i++)
3647  u[i] = 1./a[ia[i]];
3648 
3649  // PCG iterations
3650  std::vector<Real> r(n), p(n), z(n);
3651  info = pcg_ic0(n, ia, ja, a, u, x, b, r, p, z, eps, maxite, msglev);
3652 
3653  // Mat sparse_global;
3654  // MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,n,n,ia,ja,a,&sparse_global);
3655 
3656  _logfile << "Finished Solve " << n << std::endl;
3657 
3658  return info;
3659 }

References _logfile, pcg_ic0(), and pcg_par_check().

Referenced by minJ().

◆ vertex()

Real libMesh::VariationalMeshSmoother::vertex ( Array3D< Real > &  W,
Array2D< Real > &  F,
const Array2D< Real > &  R,
const std::vector< int > &  cell_in,
Real  epsilon,
Real  w,
int  nvert,
const std::vector< Real > &  K,
const Array2D< Real > &  H,
int  me,
Real  vol,
int  f,
Real &  Vmin,
int  adp,
const std::vector< Real > &  g,
Real  sigma 
)
private

Definition at line 3409 of file mesh_smoother_vsmoother.C.

3425 {
3426  // hessian, function, gradient
3427  Array2D<Real> Q(3, nvert);
3428  basisA(Q, nvert, K, H, me);
3429 
3430  std::vector<Real> a1(3), a2(3), a3(3);
3431  for (unsigned i=0; i<_dim; i++)
3432  for (int j=0; j<nvert; j++)
3433  {
3434  a1[i] += Q[i][j]*R[cell_in[j]][0];
3435  a2[i] += Q[i][j]*R[cell_in[j]][1];
3436  if (_dim == 3)
3437  a3[i] += Q[i][j]*R[cell_in[j]][2];
3438  }
3439 
3440  // account for adaptation
3441  Real G = 0.;
3442  if (adp < 2)
3443  G = g[0];
3444  else
3445  {
3446  G = 1.;
3447  for (unsigned i=0; i<_dim; i++)
3448  {
3449  a1[i] *= std::sqrt(g[0]);
3450  a2[i] *= std::sqrt(g[1]);
3451  if (_dim == 3)
3452  a3[i] *= std::sqrt(g[2]);
3453  }
3454  }
3455 
3456  Real
3457  det = 0.,
3458  tr = 0.,
3459  phit = 0.;
3460 
3461  std::vector<Real> av1(3), av2(3), av3(3);
3462 
3463  if (_dim == 2)
3464  {
3465  av1[0] = a2[1];
3466  av1[1] = -a2[0];
3467  av2[0] = -a1[1];
3468  av2[1] = a1[0];
3469  det = jac2(a1[0], a1[1], a2[0], a2[1]);
3470  for (int i=0; i<2; i++)
3471  tr += 0.5*(a1[i]*a1[i] + a2[i]*a2[i]);
3472 
3473  phit = (1-w)*tr + w*0.5*(vol/G + det*det*G/vol);
3474  }
3475 
3476  if (_dim == 3)
3477  {
3478  av1[0] = a2[1]*a3[2] - a2[2]*a3[1];
3479  av1[1] = a2[2]*a3[0] - a2[0]*a3[2];
3480  av1[2] = a2[0]*a3[1] - a2[1]*a3[0];
3481 
3482  av2[0] = a3[1]*a1[2] - a3[2]*a1[1];
3483  av2[1] = a3[2]*a1[0] - a3[0]*a1[2];
3484  av2[2] = a3[0]*a1[1] - a3[1]*a1[0];
3485 
3486  av3[0] = a1[1]*a2[2] - a1[2]*a2[1];
3487  av3[1] = a1[2]*a2[0] - a1[0]*a2[2];
3488  av3[2] = a1[0]*a2[1] - a1[1]*a2[0];
3489 
3490  det = jac3(a1[0], a1[1], a1[2],
3491  a2[0], a2[1], a2[2],
3492  a3[0], a3[1], a3[2]);
3493 
3494  for (int i=0; i<3; i++)
3495  tr += (1./3.)*(a1[i]*a1[i] + a2[i]*a2[i] + a3[i]*a3[i]);
3496 
3497  phit = (1-w)*pow(tr,1.5) + w*0.5*(vol/G + det*det*G/vol);
3498  }
3499 
3500  Real dchi = 0.5 + det*0.5/std::sqrt(epsilon*epsilon + det*det);
3501  Real chi = 0.5*(det + std::sqrt(epsilon*epsilon + det*det));
3502  Real fet = (chi != 0.) ? phit/chi : 1.e32;
3503 
3504  // Set return value reference
3505  qmin = det*G;
3506 
3507  // gradient and Hessian
3508  if (f == 0)
3509  {
3510  Array3D<Real> P(3, 3, 3);
3511  Array3D<Real> d2phi(3, 3, 3);
3512  Array2D<Real> dphi(3, 3);
3513  Array2D<Real> dfe(3, 3);
3514 
3515  if (_dim == 2)
3516  {
3517  for (int i=0; i<2; i++)
3518  {
3519  dphi[0][i] = (1-w)*a1[i] + w*G*det*av1[i]/vol;
3520  dphi[1][i] = (1-w)*a2[i] + w*G*det*av2[i]/vol;
3521  }
3522 
3523  for (int i=0; i<2; i++)
3524  for (int j=0; j<2; j++)
3525  {
3526  d2phi[0][i][j] = w*G*av1[i]*av1[j]/vol;
3527  d2phi[1][i][j] = w*G*av2[i]*av2[j]/vol;
3528 
3529  if (i == j)
3530  for (int k=0; k<2; k++)
3531  d2phi[k][i][j] += 1.-w;
3532  }
3533 
3534  for (int i=0; i<2; i++)
3535  {
3536  dfe[0][i] = (dphi[0][i] - fet*dchi*av1[i])/chi;
3537  dfe[1][i] = (dphi[1][i] - fet*dchi*av2[i])/chi;
3538  }
3539 
3540  for (int i=0; i<2; i++)
3541  for (int j=0; j<2; j++)
3542  {
3543  P[0][i][j] = (d2phi[0][i][j] - dchi*(av1[i]*dfe[0][j] + dfe[0][i]*av1[j]))/chi;
3544  P[1][i][j] = (d2phi[1][i][j] - dchi*(av2[i]*dfe[1][j] + dfe[1][i]*av2[j]))/chi;
3545  }
3546  }
3547 
3548  if (_dim == 3)
3549  {
3550  for (int i=0; i<3; i++)
3551  {
3552  dphi[0][i] = (1-w)*std::sqrt(tr)*a1[i] + w*G*det*av1[i]/vol;
3553  dphi[1][i] = (1-w)*std::sqrt(tr)*a2[i] + w*G*det*av2[i]/vol;
3554  dphi[2][i] = (1-w)*std::sqrt(tr)*a3[i] + w*G*det*av3[i]/vol;
3555  }
3556 
3557  for (int i=0; i<3; i++)
3558  {
3559  if (tr != 0)
3560  {
3561  d2phi[0][i][i] = (1-w)/(3.*std::sqrt(tr))*a1[i]*a1[i] + w*G*av1[i]*av1[i]/vol;
3562  d2phi[1][i][i] = (1-w)/(3.*std::sqrt(tr))*a2[i]*a2[i] + w*G*av2[i]*av2[i]/vol;
3563  d2phi[2][i][i] = (1-w)/(3.*std::sqrt(tr))*a3[i]*a3[i] + w*G*av3[i]*av3[i]/vol;
3564  }
3565  else
3566  {
3567  for (int k=0; k<3; k++)
3568  d2phi[k][i][i] = 0.;
3569  }
3570 
3571  for (int k=0; k<3; k++)
3572  d2phi[k][i][i] += (1-w)*std::sqrt(tr);
3573  }
3574 
3575  const Real con = 100.;
3576 
3577  for (int i=0; i<3; i++)
3578  {
3579  dfe[0][i] = (dphi[0][i] - fet*dchi*av1[i] + con*epsilon*a1[i])/chi;
3580  dfe[1][i] = (dphi[1][i] - fet*dchi*av2[i] + con*epsilon*a2[i])/chi;
3581  dfe[2][i] = (dphi[2][i] - fet*dchi*av3[i] + con*epsilon*a3[i])/chi;
3582  }
3583 
3584  for (int i=0; i<3; i++)
3585  {
3586  P[0][i][i] = (d2phi[0][i][i] - dchi*(av1[i]*dfe[0][i] + dfe[0][i]*av1[i]) + con*epsilon)/chi;
3587  P[1][i][i] = (d2phi[1][i][i] - dchi*(av2[i]*dfe[1][i] + dfe[1][i]*av2[i]) + con*epsilon)/chi;
3588  P[2][i][i] = (d2phi[2][i][i] - dchi*(av3[i]*dfe[2][i] + dfe[2][i]*av3[i]) + con*epsilon)/chi;
3589  }
3590 
3591  for (int k=0; k<3; k++)
3592  for (int i=0; i<3; i++)
3593  for (int j=0; j<3; j++)
3594  if (i != j)
3595  P[k][i][j] = 0.;
3596  }
3597 
3598  /*--------matrix W, right side F---------------------*/
3599  Array2D<Real> gpr(3, nvert);
3600 
3601  for (unsigned i1=0; i1<_dim; i1++)
3602  {
3603  for (unsigned i=0; i<_dim; i++)
3604  for (int j=0; j<nvert; j++)
3605  for (unsigned k=0; k<_dim; k++)
3606  gpr[i][j] += P[i1][i][k]*Q[k][j];
3607 
3608  for (int i=0; i<nvert; i++)
3609  for (int j=0; j<nvert; j++)
3610  for (unsigned k=0; k<_dim; k++)
3611  W[i1][i][j] += Q[k][i]*gpr[k][j]*sigma;
3612 
3613  for (int i=0; i<nvert; i++)
3614  for (int k=0; k<2; k++)
3615  F[i1][i] += Q[k][i]*dfe[i1][k]*sigma;
3616  }
3617  }
3618 
3619  return fet*sigma;
3620 }

References _dim, basisA(), jac2(), jac3(), std::pow(), libMesh::Real, and std::sqrt().

Referenced by localP().

◆ writegr()

int libMesh::VariationalMeshSmoother::writegr ( const Array2D< Real > &  R)
private

Definition at line 223 of file mesh_smoother_vsmoother.C.

224 {
225  libMesh::out << "Starting writegr" << std::endl;
226 
227  // Adjust nodal coordinates to new positions
228  {
229  libmesh_assert_equal_to(_dist_norm, 0.);
230  _dist_norm = 0;
231  int i = 0;
232  for (auto & node : _mesh.node_ptr_range())
233  {
234  Real total_dist = 0.;
235 
236  // Get a reference to the node
237  Node & node_ref = *node;
238 
239  // For each node set its X Y [Z] coordinates
240  for (unsigned int j=0; j<_dim; j++)
241  {
242  Real distance = R[i][j] - node_ref(j);
243 
244  // Save the squares of the distance
245  total_dist += Utility::pow<2>(distance);
246 
247  node_ref(j) += distance * _percent_to_move;
248  }
249 
250  libmesh_assert_greater_equal (total_dist, 0.);
251 
252  // Add the distance this node moved to the global distance
253  _dist_norm += total_dist;
254 
255  i++;
256  }
257 
258  // Relative "error"
259  _dist_norm = std::sqrt(_dist_norm/_mesh.n_nodes());
260  }
261 
262  libMesh::out << "Finished writegr" << std::endl;
263  return 0;
264 }

References _dim, _dist_norm, libMesh::MeshSmoother::_mesh, _percent_to_move, distance(), libMesh::MeshBase::n_nodes(), libMesh::MeshBase::node_ptr_range(), libMesh::out, libMesh::Real, and std::sqrt().

Referenced by smooth().

Member Data Documentation

◆ _adapt_data

std::vector<float>* libMesh::VariationalMeshSmoother::_adapt_data
private

Vector for holding adaptive data.

Definition at line 174 of file mesh_smoother_vsmoother.h.

Referenced by adapt_minimum(), adjust_adapt_data(), and read_adp().

◆ _adaptive_func

AdaptType libMesh::VariationalMeshSmoother::_adaptive_func
private

Definition at line 184 of file mesh_smoother_vsmoother.h.

Referenced by smooth().

◆ _area_of_interest

const UnstructuredMesh* libMesh::VariationalMeshSmoother::_area_of_interest
private

Area of Interest Mesh.

Definition at line 217 of file mesh_smoother_vsmoother.h.

Referenced by adjust_adapt_data(), and read_adp().

◆ _dim

const unsigned libMesh::VariationalMeshSmoother::_dim
private

Smoother control variables.

Definition at line 179 of file mesh_smoother_vsmoother.h.

Referenced by adp_renew(), avertex(), basisA(), localP(), maxE(), metr_data_gen(), minJ(), minq(), readgr(), readmetr(), smooth(), vertex(), and writegr().

◆ _dist_norm

Real libMesh::VariationalMeshSmoother::_dist_norm
private

Records a relative "distance moved".

Definition at line 164 of file mesh_smoother_vsmoother.h.

Referenced by smooth(), and writegr().

◆ _distance

Real libMesh::VariationalMeshSmoother::_distance
private

Max distance of the last set of movement.

Definition at line 154 of file mesh_smoother_vsmoother.h.

Referenced by distance_moved(), and smooth().

◆ _generate_data

bool libMesh::VariationalMeshSmoother::_generate_data
private

Definition at line 186 of file mesh_smoother_vsmoother.h.

Referenced by set_generate_data(), and smooth().

◆ _hanging_nodes

std::map<dof_id_type, std::vector<dof_id_type> > libMesh::VariationalMeshSmoother::_hanging_nodes
private

Map for hanging_nodes.

Definition at line 169 of file mesh_smoother_vsmoother.h.

Referenced by readgr(), and smooth().

◆ _logfile

std::ofstream libMesh::VariationalMeshSmoother::_logfile
private

All output (including debugging) is sent to the _logfile.

Definition at line 212 of file mesh_smoother_vsmoother.h.

Referenced by full_smooth(), minJ(), minJ_BC(), pcg_ic0(), pcg_par_check(), smooth(), and solver().

◆ _maxiter

const unsigned libMesh::VariationalMeshSmoother::_maxiter
private

Definition at line 181 of file mesh_smoother_vsmoother.h.

Referenced by smooth().

◆ _mesh

UnstructuredMesh& libMesh::MeshSmoother::_mesh
protectedinherited

Definition at line 61 of file mesh_smoother.h.

Referenced by adjust_adapt_data(), libMesh::LaplaceMeshSmoother::allgather_graph(), libMesh::LaplaceMeshSmoother::init(), metr_data_gen(), readgr(), libMesh::LaplaceMeshSmoother::smooth(), smooth(), and writegr().

◆ _metric

MetricType libMesh::VariationalMeshSmoother::_metric
private

Definition at line 183 of file mesh_smoother_vsmoother.h.

Referenced by set_metric(), and smooth().

◆ _miniter

const unsigned libMesh::VariationalMeshSmoother::_miniter
private

Definition at line 180 of file mesh_smoother_vsmoother.h.

Referenced by smooth().

◆ _miniterBC

const unsigned libMesh::VariationalMeshSmoother::_miniterBC
private

Definition at line 182 of file mesh_smoother_vsmoother.h.

Referenced by smooth().

◆ _n_cells

dof_id_type libMesh::VariationalMeshSmoother::_n_cells
private

The number of active elements in the Mesh at the time of smoothing.

Not set until smooth() is actually called to mimic the original code's behavior.

Definition at line 200 of file mesh_smoother_vsmoother.h.

Referenced by adp_renew(), full_smooth(), maxE(), metr_data_gen(), minJ(), minJ_BC(), minq(), readmetr(), and smooth().

◆ _n_hanging_edges

dof_id_type libMesh::VariationalMeshSmoother::_n_hanging_edges
private

The number of hanging node edges in the Mesh at the time of smoothing.

Not set until smooth() is actually called to mimic the original code's behavior.

Definition at line 207 of file mesh_smoother_vsmoother.h.

Referenced by full_smooth(), minJ(), and smooth().

◆ _n_nodes

dof_id_type libMesh::VariationalMeshSmoother::_n_nodes
private

The number of nodes in the Mesh at the time of smoothing.

Not set until smooth() is actually called to mimic the original code's behavior.

Definition at line 193 of file mesh_smoother_vsmoother.h.

Referenced by adp_renew(), full_smooth(), metr_data_gen(), minJ(), minJ_BC(), and smooth().

◆ _percent_to_move

const Real libMesh::VariationalMeshSmoother::_percent_to_move
private

Dampening factor.

Definition at line 159 of file mesh_smoother_vsmoother.h.

Referenced by writegr().

◆ _theta

const Real libMesh::VariationalMeshSmoother::_theta
private

Definition at line 185 of file mesh_smoother_vsmoother.h.

Referenced by smooth().

The documentation for this class was generated from the following files:
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:67
libMesh::pi
const Real pi
.
Definition: libmesh.h:237
libMesh::VariationalMeshSmoother::avertex
Real avertex(const std::vector< Real > &afun, std::vector< Real > &G, const Array2D< Real > &R, const std::vector< int > &cell_in, int nvert, int adp)
Definition: mesh_smoother_vsmoother.C:3317
libMesh::VariationalMeshSmoother::pcg_par_check
int pcg_par_check(int n, const std::vector< int > &ia, const std::vector< int > &ja, const std::vector< Real > &a, Real eps, int maxite, int msglev)
Definition: mesh_smoother_vsmoother.C:3664
libMesh::MeshSmoother::_mesh
UnstructuredMesh & _mesh
Definition: mesh_smoother.h:61
libMesh::VariationalMeshSmoother::adp_renew
void adp_renew(const Array2D< Real > &R, const Array2D< int > &cells, std::vector< Real > &afun, int adp)
Definition: mesh_smoother_vsmoother.C:829
libMesh::VariationalMeshSmoother::_n_nodes
dof_id_type _n_nodes
The number of nodes in the Mesh at the time of smoothing.
Definition: mesh_smoother_vsmoother.h:193
libMesh::VariationalMeshSmoother::_logfile
std::ofstream _logfile
All output (including debugging) is sent to the _logfile.
Definition: mesh_smoother_vsmoother.h:212
libMesh::VariationalMeshSmoother::_area_of_interest
const UnstructuredMesh * _area_of_interest
Area of Interest Mesh.
Definition: mesh_smoother_vsmoother.h:217
libMesh::VariationalMeshSmoother::smooth
virtual void smooth() override
Redefinition of the smooth function from the base class.
Definition: mesh_smoother_vsmoother.h:125
libMesh::MeshBase::active_element_ptr_range
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
libMesh::VariationalMeshSmoother::readmetr
int readmetr(std::string name, Array3D< Real > &H)
Definition: mesh_smoother_vsmoother.C:490
end
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end.
Definition: variant_filter_iterator.h:343
std::sqrt
MetaPhysicL::DualNumber< T, D > sqrt(const MetaPhysicL::DualNumber< T, D > &in)
libMesh::VariationalMeshSmoother::NONE
Definition: mesh_smoother_vsmoother.h:110
libMesh::VariationalMeshSmoother::UNIFORM
Definition: mesh_smoother_vsmoother.h:102
libMesh::VariationalMeshSmoother::readgr
int readgr(Array2D< Real > &R, std::vector< int > &mask, Array2D< int > &cells, std::vector< int > &mcells, std::vector< int > &edges, std::vector< int > &hnodes)
Definition: mesh_smoother_vsmoother.C:269
mesh
MeshBase & mesh
Definition: mesh_communication.C:1257
libMesh::MeshTools::build_nodes_to_elem_map
void build_nodes_to_elem_map(const MeshBase &mesh, std::vector< std::vector< dof_id_type >> &nodes_to_elem_map)
After calling this function the input vector nodes_to_elem_map will contain the node to element conne...
Definition: mesh_tools.C:248
libMesh::VariationalMeshSmoother::_n_hanging_edges
dof_id_type _n_hanging_edges
The number of hanging node edges in the Mesh at the time of smoothing.
Definition: mesh_smoother_vsmoother.h:207
libMesh::VariationalMeshSmoother::writegr
int writegr(const Array2D< Real > &R)
Definition: mesh_smoother_vsmoother.C:223
libMesh::VariationalMeshSmoother::minJ
Real minJ(Array2D< Real > &R, const std::vector< int > &mask, const Array2D< int > &cells, const std::vector< int > &mcells, Real epsilon, Real w, int me, const Array3D< Real > &H, Real vol, const std::vector< int > &edges, const std::vector< int > &hnodes, int msglev, Real &Vmin, Real &emax, Real &qmin, int adp, const std::vector< Real > &afun)
Definition: mesh_smoother_vsmoother.C:1887
libMesh::MeshBase::elements_begin
virtual element_iterator elements_begin()=0
Iterate over all the elements in the Mesh.
libMesh::VariationalMeshSmoother::_hanging_nodes
std::map< dof_id_type, std::vector< dof_id_type > > _hanging_nodes
Map for hanging_nodes.
Definition: mesh_smoother_vsmoother.h:169
libMesh::VariationalMeshSmoother::adapt_minimum
float adapt_minimum() const
Definition: mesh_smoother_vsmoother.C:512
libMesh::MeshSmoother::MeshSmoother
MeshSmoother(UnstructuredMesh &mesh)
Constructor.
Definition: mesh_smoother.h:46
libMesh::VariationalMeshSmoother::_adapt_data
std::vector< float > * _adapt_data
Vector for holding adaptive data.
Definition: mesh_smoother_vsmoother.h:174
libMesh::VariationalMeshSmoother::NODE
Definition: mesh_smoother_vsmoother.h:111
std::abs
MetaPhysicL::DualNumber< T, D > abs(const MetaPhysicL::DualNumber< T, D > &in)
libMesh::MeshBase::node_ptr_range
virtual SimpleRange< node_iterator > node_ptr_range()=0
libMesh::VariationalMeshSmoother::jac2
Real jac2(Real x1, Real y1, Real x2, Real y2)
Definition: mesh_smoother_vsmoother.C:610
libMesh::VariationalMeshSmoother::_generate_data
bool _generate_data
Definition: mesh_smoother_vsmoother.h:186
libMesh::VariationalMeshSmoother::vertex
Real vertex(Array3D< Real > &W, Array2D< Real > &F, const Array2D< Real > &R, const std::vector< int > &cell_in, Real epsilon, Real w, int nvert, const std::vector< Real > &K, const Array2D< Real > &H, int me, Real vol, int f, Real &Vmin, int adp, const std::vector< Real > &g, Real sigma)
Definition: mesh_smoother_vsmoother.C:3409
libMesh::VariationalMeshSmoother::jac3
Real jac3(Real x1, Real y1, Real z1, Real x2, Real y2, Real z2, Real x3, Real y3, Real z3)
Definition: mesh_smoother_vsmoother.C:595
libMesh::VariationalMeshSmoother::localP
Real localP(Array3D< Real > &W, Array2D< Real > &F, Array2D< Real > &R, const std::vector< int > &cell_in, const std::vector< int > &mask, Real epsilon, Real w, int nvert, const Array2D< Real > &H, int me, Real vol, int f, Real &Vmin, Real &qmin, int adp, const std::vector< Real > &afun, std::vector< Real > &Gloc)
Definition: mesh_smoother_vsmoother.C:2979
libMesh::VariationalMeshSmoother::_dist_norm
Real _dist_norm
Records a relative "distance moved".
Definition: mesh_smoother_vsmoother.h:164
std::pow
double pow(double a, int b)
Definition: libmesh_augment_std_namespace.h:58
libMesh::VariationalMeshSmoother::DIRECTIONAL
Definition: mesh_smoother_vsmoother.h:104
libMesh::VariationalMeshSmoother::maxE
Real maxE(Array2D< Real > &R, const Array2D< int > &cells, const std::vector< int > &mcells, int me, const Array3D< Real > &H, Real v, Real epsilon, Real w, std::vector< Real > &Gamma, Real &qmin)
Definition: mesh_smoother_vsmoother.C:1077
libMesh::VariationalMeshSmoother::minq
Real minq(const Array2D< Real > &R, const Array2D< int > &cells, const std::vector< int > &mcells, int me, const Array3D< Real > &H, Real &vol, Real &Vmin)
Definition: mesh_smoother_vsmoother.C:1484
A
static PetscErrorCode Mat * A
Definition: petscdmlibmeshimpl.C:1026
libMesh::MeshTools::find_nodal_neighbors
void find_nodal_neighbors(const MeshBase &mesh, const Node &n, const std::vector< std::vector< const Elem * >> &nodes_to_elem_map, std::vector< const Node * > &neighbors)
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:743
libMesh::VariationalMeshSmoother::CELL
Definition: mesh_smoother_vsmoother.h:109
libMesh::VariationalMeshSmoother::_theta
const Real _theta
Definition: mesh_smoother_vsmoother.h:185
libMesh::VariationalMeshSmoother::VOLUMETRIC
Definition: mesh_smoother_vsmoother.h:103
libMesh::VariationalMeshSmoother::_n_cells
dof_id_type _n_cells
The number of active elements in the Mesh at the time of smoothing.
Definition: mesh_smoother_vsmoother.h:200
libMesh::VariationalMeshSmoother::adjust_adapt_data
void adjust_adapt_data()
Definition: mesh_smoother_vsmoother.C:531
libMesh::MeshBase::n_nodes
virtual dof_id_type n_nodes() const =0
distance
Real distance(const Point &p)
Definition: subdomains_ex3.C:50
libMesh::VariationalMeshSmoother::minJ_BC
Real minJ_BC(Array2D< Real > &R, const std::vector< int > &mask, const Array2D< int > &cells, const std::vector< int > &mcells, Real epsilon, Real w, int me, const Array3D< Real > &H, Real vol, int msglev, Real &Vmin, Real &emax, Real &qmin, int adp, const std::vector< Real > &afun, int NCN)
Definition: mesh_smoother_vsmoother.C:2359
libMesh::VariationalMeshSmoother::solver
int solver(int n, const std::vector< int > &ia, const std::vector< int > &ja, const std::vector< Real > &a, std::vector< Real > &x, const std::vector< Real > &b, Real eps, int maxite, int msglev)
Definition: mesh_smoother_vsmoother.C:3627
libMesh::VariationalMeshSmoother::_adaptive_func
AdaptType _adaptive_func
Definition: mesh_smoother_vsmoother.h:184
libMesh::VariationalMeshSmoother::_dim
const unsigned _dim
Smoother control variables.
Definition: mesh_smoother_vsmoother.h:179
libMesh::VariationalMeshSmoother::_miniter
const unsigned _miniter
Definition: mesh_smoother_vsmoother.h:180
libMesh::MeshTools::find_hanging_nodes_and_parents
void find_hanging_nodes_and_parents(const MeshBase &mesh, std::map< dof_id_type, std::vector< dof_id_type >> &hanging_nodes)
Given a mesh hanging_nodes will be filled with an associative array keyed off the global id of all th...
Definition: mesh_tools.C:1078
libMesh::MeshBase::elements_end
virtual element_iterator elements_end()=0
libMesh::VariationalMeshSmoother::metr_data_gen
void metr_data_gen(std::string grid, std::string metr, int me)
Definition: mesh_smoother_vsmoother.C:3988
libMesh::VariationalMeshSmoother::full_smooth
void full_smooth(Array2D< Real > &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, Real w, const std::vector< int > &iter, int me, const Array3D< Real > &H, int adp, int gr)
Definition: mesh_smoother_vsmoother.C:876
libMesh::VariationalMeshSmoother::basisA
int basisA(Array2D< Real > &Q, int nvert, const std::vector< Real > &K, const Array2D< Real > &H, int me)
Definition: mesh_smoother_vsmoother.C:621
libMesh::VariationalMeshSmoother::_maxiter
const unsigned _maxiter
Definition: mesh_smoother_vsmoother.h:181
libMesh::VariationalMeshSmoother::_percent_to_move
const Real _percent_to_move
Dampening factor.
Definition: mesh_smoother_vsmoother.h:159
libMesh::MeshTools::find_boundary_nodes
void find_boundary_nodes(const MeshBase &mesh, std::vector< bool > &on_boundary)
Definition: mesh_tools.C:306
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:121
libMesh::VariationalMeshSmoother::_miniterBC
const unsigned _miniterBC
Definition: mesh_smoother_vsmoother.h:182
libMesh::VariationalMeshSmoother::read_adp
int read_adp(std::vector< Real > &afun)
Definition: mesh_smoother_vsmoother.C:573
libMesh::out
OStreamProxy out
libMesh::VariationalMeshSmoother::_distance
Real _distance
Max distance of the last set of movement.
Definition: mesh_smoother_vsmoother.h:154
libMesh::MeshBase::n_active_elem
virtual dof_id_type n_active_elem() const =0
int
void ErrorVector unsigned int
Definition: adjoints_ex3.C:360
libMesh::Quality::name
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
libMesh::VariationalMeshSmoother::_metric
MetricType _metric
Definition: mesh_smoother_vsmoother.h:183
libMesh::VariationalMeshSmoother::pcg_ic0
int pcg_ic0(int n, const std::vector< int > &ia, const std::vector< int > &ja, const std::vector< Real > &a, const std::vector< Real > &u, std::vector< Real > &x, const std::vector< Real > &b, std::vector< Real > &r, std::vector< Real > &p, std::vector< Real > &z, Real eps, int maxite, int msglev)
Definition: mesh_smoother_vsmoother.C:3806
Generated on Sat Jan 25 2020 12:08:57 for libMesh by   doxygen 1.8.16