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 ()=default
 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

◆ MetricType

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  :
52  _dist_norm(0.),
53  _adapt_data(nullptr),
54  _dim(mesh.mesh_dimension()),
55  _miniter(miniter),
56  _maxiter(maxiter),
57  _miniterBC(miniterBC),
60  _theta(theta),
61  _generate_data(false),
62  _n_nodes(0),
63  _n_cells(0),
65  _area_of_interest(nullptr)
66 {}
dof_id_type _n_hanging_edges
The number of hanging node edges in the Mesh at the time of smoothing.
MeshBase & mesh
dof_id_type _n_cells
The number of active elements in the Mesh at the time of smoothing.
dof_id_type _n_nodes
The number of nodes in the Mesh at the time of smoothing.
Real _dist_norm
Records a relative "distance moved".
std::vector< float > * _adapt_data
Vector for holding adaptive data.
MeshSmoother(UnstructuredMesh &mesh)
Constructor.
Definition: mesh_smoother.h:46
const Real _percent_to_move
Dampening factor.
const UnstructuredMesh * _area_of_interest
Area of Interest Mesh.
const unsigned _dim
Smoother control variables.

◆ 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  :
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),
88  _theta(theta),
89  _generate_data(false),
90  _n_nodes(0),
91  _n_cells(0),
93  _area_of_interest(nullptr)
94 {}
dof_id_type _n_hanging_edges
The number of hanging node edges in the Mesh at the time of smoothing.
MeshBase & mesh
dof_id_type _n_cells
The number of active elements in the Mesh at the time of smoothing.
dof_id_type _n_nodes
The number of nodes in the Mesh at the time of smoothing.
Real _dist_norm
Records a relative "distance moved".
std::vector< float > * _adapt_data
Vector for holding adaptive data.
MeshSmoother(UnstructuredMesh &mesh)
Constructor.
Definition: mesh_smoother.h:46
const Real _percent_to_move
Dampening factor.
const UnstructuredMesh * _area_of_interest
Area of Interest Mesh.
const unsigned _dim
Smoother control variables.

◆ 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  :
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),
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 {}
dof_id_type _n_hanging_edges
The number of hanging node edges in the Mesh at the time of smoothing.
MeshBase & mesh
dof_id_type _n_cells
The number of active elements in the Mesh at the time of smoothing.
dof_id_type _n_nodes
The number of nodes in the Mesh at the time of smoothing.
Real _dist_norm
Records a relative "distance moved".
std::vector< float > * _adapt_data
Vector for holding adaptive data.
MeshSmoother(UnstructuredMesh &mesh)
Constructor.
Definition: mesh_smoother.h:46
const Real _percent_to_move
Dampening factor.
const UnstructuredMesh * _area_of_interest
Area of Interest Mesh.
const unsigned _dim
Smoother control variables.

◆ ~VariationalMeshSmoother()

virtual libMesh::VariationalMeshSmoother::~VariationalMeshSmoother ( )
virtualdefault

Destructor.

Member Function Documentation

◆ adapt_minimum()

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

Definition at line 512 of file mesh_smoother_vsmoother.C.

References _adapt_data.

Referenced by adjust_adapt_data().

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 }
std::vector< float > * _adapt_data
Vector for holding adaptive data.

◆ adjust_adapt_data()

void libMesh::VariationalMeshSmoother::adjust_adapt_data ( )
private

Definition at line 531 of file mesh_smoother_vsmoother.C.

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

Referenced by read_adp().

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 avg = (*el)->vertex_average();
552  bool in_aoe = false;
553 
554  // See if the element's vertex average 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(avg))
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 }
UnstructuredMesh & _mesh
Definition: mesh_smoother.h:61
std::vector< float > * _adapt_data
Vector for holding adaptive data.
const UnstructuredMesh * _area_of_interest
Area of Interest Mesh.

◆ 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.

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

Referenced by full_smooth().

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  // This was an unused variable since Derek's first import of
864  // Larisa's code!? Could something have been miscopied?
865 /*
866  Real z = 0;
867  for (unsigned j=0; j<_dim; j++)
868  z += R[i][j];
869 */
870 
871  // adaptive function, node based
872  afun[i] = 5*std::sin(R[i][0]);
873  }
874  }
875 }
dof_id_type _n_cells
The number of active elements in the Mesh at the time of smoothing.
dof_id_type _n_nodes
The number of nodes in the Mesh at the time of smoothing.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
uint8_t dof_id_type
Definition: id_types.h:67
const unsigned _dim
Smoother control variables.

◆ 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 3323 of file mesh_smoother_vsmoother.C.

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

Referenced by localP().

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

◆ 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.

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

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

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 }
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
Definition: type_vector.h:53
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const unsigned _dim
Smoother control variables.

◆ 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.

References _distance.

137 { return _distance; }
Real _distance
Max distance of the last set of movement.

◆ 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 880 of file mesh_smoother_vsmoother.C.

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().

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

◆ gener()

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

Definition at line 3924 of file mesh_smoother_vsmoother.C.

References libMesh::Real.

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

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

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

◆ 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.

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

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

◆ 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 2985 of file mesh_smoother_vsmoother.C.

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

Referenced by minJ(), and minJ_BC().

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

◆ 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 1083 of file mesh_smoother_vsmoother.C.

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

Referenced by full_smooth().

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

◆ metr_data_gen()

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

Definition at line 3994 of file mesh_smoother_vsmoother.C.

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

Referenced by smooth().

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

◆ 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 1893 of file mesh_smoother_vsmoother.C.

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

Referenced by full_smooth().

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

◆ 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 2365 of file mesh_smoother_vsmoother.C.

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

Referenced by full_smooth().

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

◆ 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 1490 of file mesh_smoother_vsmoother.C.

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

Referenced by full_smooth().

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

◆ 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 3812 of file mesh_smoother_vsmoother.C.

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

Referenced by solver().

3825 {
3826  // Return value
3827  int i = 0;
3828 
3829  Real
3830  rhr = 0.,
3831  rhr0 = 0.,
3832  rhrold = 0.,
3833  rr0 = 0.;
3834 
3835  for (i=0; i<=maxite; i++)
3836  {
3837  if (i == 0)
3838  p = x;
3839 
3840  // z=Ap
3841  for (int ii=0; ii<n; ii++)
3842  z[ii] = 0.;
3843 
3844  for (int ii=0; ii<n; ii++)
3845  {
3846  z[ii] += a[ia[ii]]*p[ii];
3847 
3848  for (int j=ia[ii]+1; j<ia[ii+1]; j++)
3849  {
3850  z[ii] += a[j]*p[ja[j]-1];
3851  z[ja[j]-1] += a[j]*p[ii];
3852  }
3853  }
3854 
3855  if (i == 0)
3856  for (int k=0; k<n; k++)
3857  r[k] = b[k] - z[k]; // r(0) = b - Ax(0)
3858 
3859  if (i > 0)
3860  {
3861  Real pap = 0.;
3862  for (int k=0; k<n; k++)
3863  pap += p[k]*z[k];
3864 
3865  Real alpha = rhr/pap;
3866  for (int k=0; k<n; k++)
3867  {
3868  x[k] += p[k]*alpha;
3869  r[k] -= z[k]*alpha;
3870  }
3871  rhrold = rhr;
3872  }
3873 
3874  // z = H * r
3875  for (int ii=0; ii<n; ii++)
3876  z[ii] = r[ii]*u[ii];
3877 
3878  if (i == 0)
3879  p = z;
3880 
3881  rhr = 0.;
3882  Real rr = 0.;
3883  for (int k=0; k<n; k++)
3884  {
3885  rhr += r[k]*z[k];
3886  rr += r[k]*r[k];
3887  }
3888 
3889  if (i == 0)
3890  {
3891  rhr0 = rhr;
3892  rr0 = rr;
3893  }
3894 
3895  if (msglev > 0)
3896  _logfile << i
3897  << " ) rHr ="
3898  << std::sqrt(rhr/rhr0)
3899  << " rr ="
3900  << std::sqrt(rr/rr0)
3901  << std::endl;
3902 
3903  if (rr <= eps*eps*rr0)
3904  return i;
3905 
3906  if (i >= maxite)
3907  return i;
3908 
3909  if (i > 0)
3910  {
3911  Real beta = rhr/rhrold;
3912  for (int k=0; k<n; k++)
3913  p[k] = z[k] + p[k]*beta;
3914  }
3915  } // end for i
3916 
3917  return i;
3918 }
std::ofstream _logfile
All output (including debugging) is sent to the _logfile.
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
Definition: type_vector.h:53
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ 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 3670 of file mesh_smoother_vsmoother.C.

References _logfile.

Referenced by solver().

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

◆ read_adp()

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

Definition at line 573 of file mesh_smoother_vsmoother.C.

References _adapt_data, _area_of_interest, and adjust_adapt_data().

Referenced by full_smooth().

574 {
575  std::vector<float> & adapt_data = *_adapt_data;
576 
577  if (_area_of_interest)
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 }
std::vector< float > * _adapt_data
Vector for holding adaptive data.
const UnstructuredMesh * _area_of_interest
Area of Interest Mesh.

◆ 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.

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

Referenced by metr_data_gen(), and smooth().

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 =
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 : make_range(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 (within 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 : make_range(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 : make_range(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 : make_range(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 }
std::map< dof_id_type, std::vector< dof_id_type > > _hanging_nodes
Map for hanging_nodes.
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:888
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:448
UnstructuredMesh & _mesh
Definition: mesh_smoother.h:61
std::unordered_set< dof_id_type > find_boundary_nodes(const MeshBase &mesh)
Calling this function on a 2D mesh will convert all the elements to triangles.
Definition: mesh_tools.C:516
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
Definition: type_vector.h:57
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
OStreamProxy out
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:134
const unsigned _dim
Smoother control variables.
const Real pi
.
Definition: libmesh.h:274

◆ readmetr()

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

Definition at line 490 of file mesh_smoother_vsmoother.C.

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

Referenced by smooth().

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 }
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
dof_id_type _n_cells
The number of active elements in the Mesh at the time of smoothing.
uint8_t dof_id_type
Definition: id_types.h:67
const unsigned _dim
Smoother control variables.

◆ 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.

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.

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.

References _distance, and smooth().

Referenced by smooth().

125 { _distance = this->smooth(1); }
Real _distance
Max distance of the last set of movement.
virtual void smooth() override
Redefinition of the smooth function from the base class.

◆ 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.

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().

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
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 }
std::map< dof_id_type, std::vector< dof_id_type > > _hanging_nodes
Map for hanging_nodes.
virtual dof_id_type n_active_elem() const =0
dof_id_type _n_hanging_edges
The number of hanging node edges in the Mesh at the time of smoothing.
std::ofstream _logfile
All output (including debugging) is sent to the _logfile.
UnstructuredMesh & _mesh
Definition: mesh_smoother.h:61
dof_id_type _n_cells
The number of active elements in the Mesh at the time of smoothing.
dof_id_type _n_nodes
The number of nodes in the Mesh at the time of smoothing.
int readgr(Array2D< Real > &R, std::vector< int > &mask, Array2D< int > &cells, std::vector< int > &mcells, std::vector< int > &edges, std::vector< int > &hnodes)
Real _dist_norm
Records a relative "distance moved".
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:911
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)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
int readmetr(std::string name, Array3D< Real > &H)
int writegr(const Array2D< Real > &R)
virtual dof_id_type n_nodes() const =0
const unsigned _dim
Smoother control variables.
void metr_data_gen(std::string grid, std::string metr, int me)

◆ 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 3633 of file mesh_smoother_vsmoother.C.

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

Referenced by minJ().

3642 {
3643  _logfile << "Beginning Solve " << n << std::endl;
3644 
3645  // Check parameters
3646  int info = pcg_par_check(n, ia, ja, a, eps, maxite, msglev);
3647  if (info != 0)
3648  return info;
3649 
3650  // PJ preconditioner construction
3651  std::vector<Real> u(n);
3652  for (int i=0; i<n; i++)
3653  u[i] = 1./a[ia[i]];
3654 
3655  // PCG iterations
3656  std::vector<Real> r(n), p(n), z(n);
3657  info = pcg_ic0(n, ia, ja, a, u, x, b, r, p, z, eps, maxite, msglev);
3658 
3659  // Mat sparse_global;
3660  // MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,n,n,ia,ja,a,&sparse_global);
3661 
3662  _logfile << "Finished Solve " << n << std::endl;
3663 
3664  return info;
3665 }
MPI_Info info
std::ofstream _logfile
All output (including debugging) is sent to the _logfile.
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)

◆ 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 3415 of file mesh_smoother_vsmoother.C.

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

Referenced by localP().

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

◆ writegr()

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

Definition at line 223 of file mesh_smoother_vsmoother.C.

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

Referenced by smooth().

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"
260  }
261 
262  libMesh::out << "Finished writegr" << std::endl;
263  return 0;
264 }
UnstructuredMesh & _mesh
Definition: mesh_smoother.h:61
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
Definition: type_vector.h:53
Real distance(const Point &p)
Real _dist_norm
Records a relative "distance moved".
const Real _percent_to_move
Dampening factor.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
OStreamProxy out
virtual dof_id_type n_nodes() const =0
const unsigned _dim
Smoother control variables.

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

◆ _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: