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

#include <mesh_smoother_vsmoother.h>

Inheritance diagram for libMesh::VariationalMeshSmoother:
Inheritance graph
[legend]

Classes

struct  Array2D
 
struct  Array3D
 

Public Types

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

Public Member Functions

 VariationalMeshSmoother (UnstructuredMesh &mesh, double theta=0.5, unsigned miniter=2, unsigned maxiter=5, unsigned miniterBC=5)
 
 VariationalMeshSmoother (UnstructuredMesh &mesh, std::vector< float > *adapt_data, double theta=0.5, unsigned miniter=2, unsigned maxiter=5, unsigned miniterBC=5, double percent_to_move=1)
 
 VariationalMeshSmoother (UnstructuredMesh &mesh, const UnstructuredMesh *area_of_interest, std::vector< float > *adapt_data, double theta=0.5, unsigned miniter=2, unsigned maxiter=5, unsigned miniterBC=5, double percent_to_move=1)
 
virtual ~VariationalMeshSmoother ()
 
virtual void smooth () override
 
double smooth (unsigned int n_iterations)
 
double distance_moved () const
 
void set_generate_data (bool b)
 
void set_metric (MetricType t)
 

Protected Attributes

UnstructuredMesh & _mesh
 

Private Member Functions

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

Private Attributes

double _distance
 
const double _percent_to_move
 
double _dist_norm
 
std::map< dof_id_type, std::vector< dof_id_type > > _hanging_nodes
 
std::vector< float > * _adapt_data
 
const unsigned _dim
 
const unsigned _miniter
 
const unsigned _maxiter
 
const unsigned _miniterBC
 
MetricType _metric
 
AdaptType _adaptive_func
 
const double _theta
 
bool _generate_data
 
dof_id_type _n_nodes
 
dof_id_type _n_cells
 
dof_id_type _n_hanging_edges
 
std::ofstream _logfile
 
const UnstructuredMesh * _area_of_interest
 

Detailed Description

This is an implementation of Larisa Branets' smoothing algorithms. The initial implementation was done by her, the adaptation to libmesh was completed by Derek Gaston. The code was heavily refactored into something more closely resembling C++ by John Peterson in 2014.

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

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

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

Author
Derek R. Gaston
Date
2006

Definition at line 63 of file mesh_smoother_vsmoother.h.

Member Enumeration Documentation

◆ AdaptType

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

Definition at line 107 of file mesh_smoother_vsmoother.h.

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

◆ MetricType

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

Definition at line 100 of file mesh_smoother_vsmoother.h.

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

Constructor & Destructor Documentation

◆ VariationalMeshSmoother() [1/3]

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

Simple constructor to use for smoothing purposes

Definition at line 45 of file mesh_smoother_vsmoother.C.

49  :
50  MeshSmoother(mesh),
51  _percent_to_move(1),
52  _dist_norm(0.),
53  _adapt_data(nullptr),
54  _dim(mesh.mesh_dimension()),
55  _miniter(miniter),
56  _maxiter(maxiter),
57  _miniterBC(miniterBC),
58  _metric(UNIFORM),
59  _adaptive_func(NONE),
60  _theta(theta),
61  _generate_data(false),
62  _n_nodes(0),
63  _n_cells(0),
64  _n_hanging_edges(0),
65  _area_of_interest(nullptr)
66 {}
mesh
MeshBase & mesh
Definition: mesh_communication.C:1245
libMesh::MeshSmoother::MeshSmoother
MeshSmoother(UnstructuredMesh &mesh)
Definition: mesh_smoother.h:46
libMesh::VariationalMeshSmoother::_area_of_interest
const UnstructuredMesh * _area_of_interest
Definition: mesh_smoother_vsmoother.h:217

◆ VariationalMeshSmoother() [2/3]

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

Slightly more complicated constructor for mesh redistribution based on adapt_data

Definition at line 71 of file mesh_smoother_vsmoother.C.

77  :
78  MeshSmoother(mesh),
79  _percent_to_move(percent_to_move),
80  _dist_norm(0.),
81  _adapt_data(adapt_data),
82  _dim(mesh.mesh_dimension()),
83  _miniter(miniter),
84  _maxiter(maxiter),
85  _miniterBC(miniterBC),
86  _metric(UNIFORM),
87  _adaptive_func(CELL),
88  _theta(theta),
89  _generate_data(false),
90  _n_nodes(0),
91  _n_cells(0),
92  _n_hanging_edges(0),
93  _area_of_interest(nullptr)
94 {}
mesh
MeshBase & mesh
Definition: mesh_communication.C:1245
libMesh::MeshSmoother::MeshSmoother
MeshSmoother(UnstructuredMesh &mesh)
Definition: mesh_smoother.h:46
libMesh::VariationalMeshSmoother::_area_of_interest
const UnstructuredMesh * _area_of_interest
Definition: mesh_smoother_vsmoother.h:217

◆ VariationalMeshSmoother() [3/3]

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

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

Definition at line 98 of file mesh_smoother_vsmoother.C.

105  :
106  MeshSmoother(mesh),
107  _percent_to_move(percent_to_move),
108  _dist_norm(0.),
109  _adapt_data(adapt_data),
110  _dim(mesh.mesh_dimension()),
111  _miniter(miniter),
112  _maxiter(maxiter),
113  _miniterBC(miniterBC),
114  _metric(UNIFORM),
115  _adaptive_func(CELL),
116  _theta(theta),
117  _generate_data(false),
118  _n_nodes(0),
119  _n_cells(0),
120  _n_hanging_edges(0),
121  _area_of_interest(area_of_interest)
122 {}
mesh
MeshBase & mesh
Definition: mesh_communication.C:1245
libMesh::MeshSmoother::MeshSmoother
MeshSmoother(UnstructuredMesh &mesh)
Definition: mesh_smoother.h:46
libMesh::VariationalMeshSmoother::_area_of_interest
const UnstructuredMesh * _area_of_interest
Definition: mesh_smoother_vsmoother.h:217

◆ ~VariationalMeshSmoother()

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

Destructor.

Definition at line 117 of file mesh_smoother_vsmoother.h.

117 {}

Member Function Documentation

◆ adapt_minimum()

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

Definition at line 512 of file mesh_smoother_vsmoother.C.

References _adapt_data, std::max(), and std::min().

Referenced by adjust_adapt_data().

513 {
514  float min = std::numeric_limits<float>::max();
515 
516  for (std::size_t i=0; i<_adapt_data->size(); i++)
517  {
518  // Only positive (or zero) values in the error vector
519  libmesh_assert_greater_equal ((*_adapt_data)[i], 0.);
520  min = std::min (min, (*_adapt_data)[i]);
521  }
522 
523  // ErrorVectors are for positive values
524  libmesh_assert_greater_equal (min, 0.);
525 
526  return min;
527 }
std::max
long double max(long double a, double b)
Definition: libmesh_augment_std_namespace.h:27
std::min
long double min(long double a, double b)
Definition: libmesh_augment_std_namespace.h:29

◆ 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, adapt_minimum(), libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), and std::min().

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 centroid = (*el)->centroid();
552  bool in_aoe = false;
553 
554  // See if the elements centroid lies in the aoe mesh
555  // for (aoe_el=aoe_mesh.elements_begin(); aoe_el != aoe_end_el; ++aoe_el)
556  // {
557  if ((*aoe_el)->contains_point(centroid))
558  in_aoe = true;
559  // }
560 
561  // If the element is not in the area of interest... then set
562  // it's adaptivity value to the minimum.
563  if (!in_aoe)
564  adapt_data[i] = min;
565  }
566  i++;
567  }
568 }
libMesh::MeshSmoother::_mesh
UnstructuredMesh & _mesh
Definition: mesh_smoother.h:61
libMesh::MeshBase::elements_begin
virtual element_iterator elements_begin()=0
libMesh::MeshBase::elements_end
virtual element_iterator elements_end()=0
libMesh::VariationalMeshSmoother::_area_of_interest
const UnstructuredMesh * _area_of_interest
Definition: mesh_smoother_vsmoother.h:217
std::min
long double min(long double a, double b)
Definition: libmesh_augment_std_namespace.h:29

◆ adp_renew()

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

Definition at line 829 of file mesh_smoother_vsmoother.C.

References _dim, _n_cells, and _n_nodes.

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  double
840  x = 0.,
841  y = 0.,
842  z = 0.;
843  int nvert = 0;
844 
845  while (cells[i][nvert] >= 0)
846  {
847  x += R[cells[i][nvert]][0];
848  y += R[cells[i][nvert]][1];
849  if (_dim == 3)
850  z += R[cells[i][nvert]][2];
851  nvert++;
852  }
853 
854  // adaptive function, cell based
855  afun[i] = 5.*(x+y+z);
856  }
857  }
858 
859  else if (adp > 0)
860  {
861  for (dof_id_type i=0; i<_n_nodes; i++)
862  {
863  double z = 0;
864  for (unsigned j=0; j<_dim; j++)
865  z += R[i][j];
866 
867  // adaptive function, node based
868  afun[i] = 5*std::sin(R[i][0]);
869  }
870  }
871 }
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:64

◆ avertex()

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

Definition at line 3317 of file mesh_smoother_vsmoother.C.

References _dim, basisA(), jac2(), and jac3().

Referenced by localP().

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

◆ basisA()

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

Definition at line 621 of file mesh_smoother_vsmoother.C.

References _dim.

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

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

◆ distance_moved()

double 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; }

◆ full_smooth()

void libMesh::VariationalMeshSmoother::full_smooth ( Array2D< double > &  R,
const std::vector< int > &  mask,
const Array2D< int > &  cells,
const std::vector< int > &  mcells,
const std::vector< int > &  edges,
const std::vector< int > &  hnodes,
double  w,
const std::vector< int > &  iter,
int  me,
const Array3D< double > &  H,
int  adp,
int  gr 
)
private

Definition at line 876 of file mesh_smoother_vsmoother.C.

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

Referenced by smooth().

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

◆ gener()

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

Definition at line 3918 of file mesh_smoother_vsmoother.C.

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

◆ jac2()

double libMesh::VariationalMeshSmoother::jac2 ( double  x1,
double  y1,
double  x2,
double  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()

double libMesh::VariationalMeshSmoother::jac3 ( double  x1,
double  y1,
double  z1,
double  x2,
double  y2,
double  z2,
double  x3,
double  y3,
double  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()

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

Definition at line 2979 of file mesh_smoother_vsmoother.C.

References _dim, avertex(), and vertex().

Referenced by minJ(), and minJ_BC().

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

◆ maxE()

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

Definition at line 1077 of file mesh_smoother_vsmoother.C.

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

Referenced by full_smooth().

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

◆ metr_data_gen()

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

Definition at line 3988 of file mesh_smoother_vsmoother.C.

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

Referenced by smooth().

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

◆ minJ()

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

Definition at line 1887 of file mesh_smoother_vsmoother.C.

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

Referenced by full_smooth().

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

◆ minJ_BC()

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

Definition at line 2359 of file mesh_smoother_vsmoother.C.

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

Referenced by full_smooth().

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

◆ minq()

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

Definition at line 1484 of file mesh_smoother_vsmoother.C.

References _dim, _n_cells, basisA(), jac2(), and jac3().

Referenced by full_smooth().

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

◆ pcg_ic0()

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

Definition at line 3806 of file mesh_smoother_vsmoother.C.

References _logfile.

Referenced by solver().

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

◆ pcg_par_check()

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

Definition at line 3664 of file mesh_smoother_vsmoother.C.

References _logfile.

Referenced by solver().

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

◆ read_adp()

int libMesh::VariationalMeshSmoother::read_adp ( std::vector< double > &  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)
578  adjust_adapt_data();
579 
580  std::size_t m = adapt_data.size();
581 
582  std::size_t j =0 ;
583  for (std::size_t i=0; i<m; i++)
584  if (adapt_data[i] != 0)
585  {
586  afun[j] = adapt_data[i];
587  j++;
588  }
589 
590  return 0;
591 }
libMesh::VariationalMeshSmoother::_area_of_interest
const UnstructuredMesh * _area_of_interest
Definition: mesh_smoother_vsmoother.h:217

◆ readgr()

int libMesh::VariationalMeshSmoother::readgr ( Array2D< double > &  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::MeshBase::active_element_ptr_range(), libMesh::MeshTools::build_nodes_to_elem_map(), end, libMesh::MeshTools::find_boundary_nodes(), libMesh::MeshTools::find_nodal_neighbors(), libMesh::MeshBase::node_ptr_range(), libMesh::out, libMesh::pi, and libMesh::Real.

Referenced by metr_data_gen(), and smooth().

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

◆ readmetr()

int libMesh::VariationalMeshSmoother::readmetr ( std::string  name,
Array3D< double > &  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 }
libMesh::Quality::name
std::string name(const ElemQuality q)
Definition: elem_quality.C:42
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:64

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

142 { _generate_data = b; }

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

147 { _metric = t; }

◆ 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); }

◆ smooth() [2/2]

double 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(), 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  double theta = _theta;
143 
144  // Metric file name
145  std::string metric_filename = "smoother.metric";
146  if (gr == 0 && me > 1)
147  {
148  // grid filename
149  std::string grid_filename = "smoother.grid";
150 
151  // generate metric from initial mesh (me = 2,3)
152  metr_data_gen(grid_filename, metric_filename, me);
153  }
154 
155  // Initialize the _n_nodes and _n_cells member variables
156  this->_n_nodes = _mesh.n_nodes();
157  this->_n_cells = _mesh.n_active_elem();
158 
159  // Initialize the _n_hanging_edges member variable
160  MeshTools::find_hanging_nodes_and_parents(_mesh, _hanging_nodes);
161  this->_n_hanging_edges =
162  cast_int<dof_id_type>(_hanging_nodes.size());
163 
164  std::vector<int>
165  mask(_n_nodes),
166  edges(2*_n_hanging_edges),
167  mcells(_n_cells),
168  hnodes(_n_hanging_edges);
169 
170  Array2D<double> R(_n_nodes, _dim);
171  Array2D<int> cells(_n_cells, 3*_dim + _dim%2);
172  Array3D<double> H(_n_cells, _dim, _dim);
173 
174  // initial grid
175  int vms_err = readgr(R, mask, cells, mcells, edges, hnodes);
176  if (vms_err < 0)
177  {
178  _logfile << "Error reading input mesh file" << std::endl;
179  return _dist_norm;
180  }
181 
182  if (me > 1)
183  vms_err = readmetr(metric_filename, H);
184 
185  if (vms_err < 0)
186  {
187  _logfile << "Error reading metric file" << std::endl;
188  return _dist_norm;
189  }
190 
191  std::vector<int> iter(4);
192  iter[0] = miniter;
193  iter[1] = maxiter;
194  iter[2] = miniterBC;
195 
196  // grid optimization
197  _logfile << "Starting Grid Optimization" << std::endl;
198  clock_t ticks1 = clock();
199  full_smooth(R, mask, cells, mcells, edges, hnodes, theta, iter, me, H, adp, gr);
200  clock_t ticks2 = clock();
201  _logfile << "full_smooth took ("
202  << ticks2
203  << "-"
204  << ticks1
205  << ")/"
206  << CLOCKS_PER_SEC
207  << " = "
208  << static_cast<double>(ticks2-ticks1)/static_cast<double>(CLOCKS_PER_SEC)
209  << " seconds"
210  << std::endl;
211 
212  // save result
213  _logfile << "Saving Result" << std::endl;
214  writegr(R);
215 
216  libmesh_assert_greater (_dist_norm, 0);
217  return _dist_norm;
218 }
libMesh::VariationalMeshSmoother::_hanging_nodes
std::map< dof_id_type, std::vector< dof_id_type > > _hanging_nodes
Definition: mesh_smoother_vsmoother.h:169
libMesh::VariationalMeshSmoother::full_smooth
void full_smooth(Array2D< double > &R, const std::vector< int > &mask, const Array2D< int > &cells, const std::vector< int > &mcells, const std::vector< int > &edges, const std::vector< int > &hnodes, double w, const std::vector< int > &iter, int me, const Array3D< double > &H, int adp, int gr)
Definition: mesh_smoother_vsmoother.C:876
libMesh::MeshBase::n_active_elem
virtual dof_id_type n_active_elem() const =0
libMesh::MeshSmoother::_mesh
UnstructuredMesh & _mesh
Definition: mesh_smoother.h:61
libMesh::VariationalMeshSmoother::readgr
int readgr(Array2D< double > &R, std::vector< int > &mask, Array2D< int > &cells, std::vector< int > &mcells, std::vector< int > &edges, std::vector< int > &hnodes)
Definition: mesh_smoother_vsmoother.C:269
libMesh::VariationalMeshSmoother::readmetr
int readmetr(std::string name, Array3D< double > &H)
Definition: mesh_smoother_vsmoother.C:490
libMesh::MeshTools::find_hanging_nodes_and_parents
void find_hanging_nodes_and_parents(const MeshBase &mesh, std::map< dof_id_type, std::vector< dof_id_type >> &hanging_nodes)
Definition: mesh_tools.C:1083
libMesh::VariationalMeshSmoother::writegr
int writegr(const Array2D< double > &R)
Definition: mesh_smoother_vsmoother.C:223
libMesh::MeshBase::n_nodes
virtual dof_id_type n_nodes() const =0
libMesh::VariationalMeshSmoother::metr_data_gen
void metr_data_gen(std::string grid, std::string metr, int me)
Definition: mesh_smoother_vsmoother.C:3988

◆ solver()

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

Definition at line 3627 of file mesh_smoother_vsmoother.C.

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

Referenced by minJ().

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

◆ vertex()

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

Definition at line 3409 of file mesh_smoother_vsmoother.C.

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

Referenced by localP().

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

◆ writegr()

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

Definition at line 223 of file mesh_smoother_vsmoother.C.

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

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  double total_dist = 0.;
235 
236  // Get a reference to the node
237  Node & node_ref = *node;
238 
239  // For each node set its X Y [Z] coordinates
240  for (unsigned int j=0; j<_dim; j++)
241  {
242  double distance = R[i][j] - node_ref(j);
243 
244  // Save the squares of the distance
245  total_dist += Utility::pow<2>(distance);
246 
247  node_ref(j) += distance * _percent_to_move;
248  }
249 
250  libmesh_assert_greater_equal (total_dist, 0.);
251 
252  // Add the distance this node moved to the global distance
253  _dist_norm += total_dist;
254 
255  i++;
256  }
257 
258  // Relative "error"
259  _dist_norm = std::sqrt(_dist_norm/_mesh.n_nodes());
260  }
261 
262  libMesh::out << "Finished writegr" << std::endl;
263  return 0;
264 }
libMesh::MeshSmoother::_mesh
UnstructuredMesh & _mesh
Definition: mesh_smoother.h:61
libMesh::MeshBase::node_ptr_range
virtual SimpleRange< node_iterator > node_ptr_range()=0
libMesh::out
OStreamProxy out(std::cout)
Definition: libmesh_common.h:228
libMesh::MeshBase::n_nodes
virtual dof_id_type n_nodes() const =0

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

double 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

double libMesh::VariationalMeshSmoother::_distance
private

Max distance of the last set of movement.

Definition at line 154 of file mesh_smoother_vsmoother.h.

Referenced by distance_moved(), and smooth().

◆ _generate_data

bool libMesh::VariationalMeshSmoother::_generate_data
private

Definition at line 186 of file mesh_smoother_vsmoother.h.

Referenced by set_generate_data(), and smooth().

◆ _hanging_nodes

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

Map for hanging_nodes

Definition at line 169 of file mesh_smoother_vsmoother.h.

Referenced by readgr(), and smooth().

◆ _logfile

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

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

Definition at line 212 of file mesh_smoother_vsmoother.h.

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

◆ _maxiter

const unsigned libMesh::VariationalMeshSmoother::_maxiter
private

Definition at line 181 of file mesh_smoother_vsmoother.h.

Referenced by smooth().

◆ _mesh

UnstructuredMesh& libMesh::MeshSmoother::_mesh
protectedinherited

Definition at line 61 of file mesh_smoother.h.

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

◆ _metric

MetricType libMesh::VariationalMeshSmoother::_metric
private

Definition at line 183 of file mesh_smoother_vsmoother.h.

Referenced by set_metric(), and smooth().

◆ _miniter

const unsigned libMesh::VariationalMeshSmoother::_miniter
private

Definition at line 180 of file mesh_smoother_vsmoother.h.

Referenced by smooth().

◆ _miniterBC

const unsigned libMesh::VariationalMeshSmoother::_miniterBC
private

Definition at line 182 of file mesh_smoother_vsmoother.h.

Referenced by smooth().

◆ _n_cells

dof_id_type libMesh::VariationalMeshSmoother::_n_cells
private

The number of active elements in the Mesh at the time of smoothing. Not set until smooth() is actually called to mimic the original code's behavior.

Definition at line 200 of file mesh_smoother_vsmoother.h.

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

◆ _n_hanging_edges

dof_id_type libMesh::VariationalMeshSmoother::_n_hanging_edges
private

The number of hanging node edges in the Mesh at the time of smoothing. Not set until smooth() is actually called to mimic the original code's behavior.

Definition at line 207 of file mesh_smoother_vsmoother.h.

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

◆ _n_nodes

dof_id_type libMesh::VariationalMeshSmoother::_n_nodes
private

The number of nodes in the Mesh at the time of smoothing. Not set until smooth() is actually called to mimic the original code's behavior.

Definition at line 193 of file mesh_smoother_vsmoother.h.

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

◆ _percent_to_move

const double libMesh::VariationalMeshSmoother::_percent_to_move
private

Dampening factor

Definition at line 159 of file mesh_smoother_vsmoother.h.

Referenced by writegr().

◆ _theta

const double 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:
generated by doxygen