libMesh::VariationalMeshSmoother Class Reference

#include <mesh_smoother_vsmoother.h>

Inheritance diagram for libMesh::VariationalMeshSmoother:

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

Constructor & Destructor Documentation

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

48  :
51  _dist_norm(0.),
53  _dim(mesh.mesh_dimension()),
54  _miniter(miniter),
55  _maxiter(maxiter),
56  _miniterBC(miniterBC),
59  _theta(theta),
60  _generate_data(false),
61  _n_nodes(0),
62  _n_cells(0),
65 {}
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
MeshSmoother(UnstructuredMesh &mesh)
Definition: mesh_smoother.h:46
const UnstructuredMesh * _area_of_interest
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 70 of file mesh_smoother_vsmoother.C.

76  :
78  _percent_to_move(percent_to_move),
79  _dist_norm(0.),
80  _adapt_data(adapt_data),
81  _dim(mesh.mesh_dimension()),
82  _miniter(miniter),
83  _maxiter(maxiter),
84  _miniterBC(miniterBC),
87  _theta(theta),
88  _generate_data(false),
89  _n_nodes(0),
90  _n_cells(0),
93 {}
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
MeshSmoother(UnstructuredMesh &mesh)
Definition: mesh_smoother.h:46
const UnstructuredMesh * _area_of_interest
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 97 of file mesh_smoother_vsmoother.C.

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

Destructor.

Definition at line 117 of file mesh_smoother_vsmoother.h.

117 {}

Member Function Documentation

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

Definition at line 524 of file mesh_smoother_vsmoother.C.

References _adapt_data, and std::min().

Referenced by adjust_adapt_data().

525 {
526  float min = 1.e30;
527 
528  for (std::size_t i=0; i<_adapt_data->size(); i++)
529  {
530  // Only positive (or zero) values in the error vector
531  libmesh_assert_greater_equal ((*_adapt_data)[i], 0.);
532  min = std::min (min, (*_adapt_data)[i]);
533  }
534 
535  // ErrorVectors are for positive values
536  libmesh_assert_greater_equal (min, 0.);
537 
538  return min;
539 }
long double min(long double a, double b)
void libMesh::VariationalMeshSmoother::adjust_adapt_data ( )
private

Definition at line 543 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().

544 {
545  // For convenience
546  const UnstructuredMesh & aoe_mesh = *_area_of_interest;
547  std::vector<float> & adapt_data = *_adapt_data;
548 
549  float min = adapt_minimum();
550 
551  MeshBase::const_element_iterator el = _mesh.elements_begin();
552  const MeshBase::const_element_iterator end_el = _mesh.elements_end();
553 
554  MeshBase::const_element_iterator aoe_el = aoe_mesh.elements_begin();
555  const MeshBase::const_element_iterator aoe_end_el = aoe_mesh.elements_end();
556 
557  // Counter to keep track of which spot in adapt_data we're looking at
558  for (int i=0; el!=end_el; el++)
559  {
560  // Only do this for active elements
561  if (adapt_data[i])
562  {
563  Point centroid = (*el)->centroid();
564  bool in_aoe = false;
565 
566  // See if the elements centroid lies in the aoe mesh
567  // for (aoe_el=aoe_mesh.elements_begin(); aoe_el != aoe_end_el; ++aoe_el)
568  // {
569  if ((*aoe_el)->contains_point(centroid))
570  in_aoe = true;
571  // }
572 
573  // If the element is not in the area of interest... then set
574  // it's adaptivity value to the minimum.
575  if (!in_aoe)
576  adapt_data[i] = min;
577  }
578  i++;
579  }
580 }
UnstructuredMesh & _mesh
Definition: mesh_smoother.h:61
virtual element_iterator elements_begin()=0
virtual element_iterator elements_end()=0
const UnstructuredMesh * _area_of_interest
long double min(long double a, double b)
void libMesh::VariationalMeshSmoother::adp_renew ( const Array2D< double > &  R,
const Array2D< int > &  cells,
std::vector< double > &  afun,
int  adp 
)
private

Definition at line 841 of file mesh_smoother_vsmoother.C.

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

Referenced by full_smooth().

845 {
846  // evaluates given adaptive function on the provided mesh
847  if (adp < 0)
848  {
849  for (dof_id_type i=0; i<_n_cells; i++)
850  {
851  double
852  x = 0.,
853  y = 0.,
854  z = 0.;
855  int nvert = 0;
856 
857  while (cells[i][nvert] >= 0)
858  {
859  x += R[cells[i][nvert]][0];
860  y += R[cells[i][nvert]][1];
861  if (_dim == 3)
862  z += R[cells[i][nvert]][2];
863  nvert++;
864  }
865 
866  // adaptive function, cell based
867  afun[i] = 5.*(x+y+z);
868  }
869  }
870 
871  else if (adp > 0)
872  {
873  for (dof_id_type i=0; i<_n_nodes; i++)
874  {
875  double z = 0;
876  for (unsigned j=0; j<_dim; j++)
877  z += R[i][j];
878 
879  // adaptive function, node based
880  afun[i] = 5*std::sin(R[i][0]);
881  }
882  }
883 }
PetscErrorCode Vec x
uint8_t dof_id_type
Definition: id_types.h:64
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 3309 of file mesh_smoother_vsmoother.C.

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

Referenced by localP().

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

Definition at line 633 of file mesh_smoother_vsmoother.C.

References _dim.

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

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

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

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

Definition at line 3910 of file mesh_smoother_vsmoother.C.

References libMesh::x.

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

Definition at line 622 of file mesh_smoother_vsmoother.C.

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

626 {
627  return x1*y2 - x2*y1;
628 }
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 607 of file mesh_smoother_vsmoother.C.

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

616 {
617  return x1*(y2*z3 - y3*z2) + y1*(z2*x3 - z3*x2) + z1*(x2*y3 - x3*y2);
618 }
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 2978 of file mesh_smoother_vsmoother.C.

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

Referenced by minJ(), and minJ_BC().

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

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

Referenced by full_smooth().

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

Definition at line 3980 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().

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

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

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

Referenced by full_smooth().

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

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

Referenced by full_smooth().

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

References _logfile, and libMesh::x.

Referenced by solver().

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

References _logfile.

Referenced by solver().

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

Definition at line 585 of file mesh_smoother_vsmoother.C.

References _adapt_data, _area_of_interest, and adjust_adapt_data().

Referenced by full_smooth().

586 {
587  std::vector<float> & adapt_data = *_adapt_data;
588 
589  if (_area_of_interest)
591 
592  std::size_t m = adapt_data.size();
593 
594  std::size_t j =0 ;
595  for (std::size_t i=0; i<m; i++)
596  if (adapt_data[i] != 0)
597  {
598  afun[j] = adapt_data[i];
599  j++;
600  }
601 
602  return 0;
603 }
const UnstructuredMesh * _area_of_interest
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 270 of file mesh_smoother_vsmoother.C.

References _dim, _hanging_nodes, libMesh::MeshSmoother::_mesh, std::abs(), libMesh::MeshBase::active_elements_begin(), libMesh::MeshBase::active_elements_end(), libMesh::MeshTools::build_nodes_to_elem_map(), end, libMesh::MeshTools::find_boundary_nodes(), libMesh::MeshTools::find_nodal_neighbors(), libMesh::Elem::n_vertices(), libMesh::Elem::node_id(), libMesh::MeshBase::nodes_begin(), libMesh::MeshBase::nodes_end(), libMesh::out, libMesh::pi, libMesh::Real, and libMesh::x.

Referenced by metr_data_gen(), and smooth().

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

Definition at line 502 of file mesh_smoother_vsmoother.C.

References _dim, and _n_cells.

Referenced by smooth().

504 {
505  std::ifstream infile(name.c_str());
506  std::string dummy;
507 
508  for (dof_id_type i=0; i<_n_cells; i++)
509  for (unsigned j=0; j<_dim; j++)
510  {
511  for (unsigned k=0; k<_dim; k++)
512  infile >> H[i][j][k];
513 
514  // Read to end of line and discard
515  std::getline(infile, dummy);
516  }
517 
518  return 0;
519 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:39
uint8_t dof_id_type
Definition: id_types.h:64
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.

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.

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

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

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

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

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

Referenced by minJ().

3628 {
3629  _logfile << "Beginning Solve " << n << std::endl;
3630 
3631  // Check parameters
3632  int info = pcg_par_check(n, ia, ja, a, eps, maxite, msglev);
3633  if (info != 0)
3634  return info;
3635 
3636  // PJ preconditioner construction
3637  std::vector<double> u(n);
3638  for (int i=0; i<n; i++)
3639  u[i] = 1./a[ia[i]];
3640 
3641  // PCG iterations
3642  std::vector<double> r(n), p(n), z(n);
3643  info = pcg_ic0(n, ia, ja, a, u, x, b, r, p, z, eps, maxite, msglev);
3644 
3645  // Mat sparse_global;
3646  // MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,n,n,ia,ja,a,&sparse_global);
3647 
3648  _logfile << "Finished Solve " << n << std::endl;
3649 
3650  return info;
3651 }
PetscErrorCode Vec x
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)
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)
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 3401 of file mesh_smoother_vsmoother.C.

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

Referenced by localP().

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

Definition at line 222 of file mesh_smoother_vsmoother.C.

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

Referenced by smooth().

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

Member Data Documentation

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

AdaptType libMesh::VariationalMeshSmoother::_adaptive_func
private

Definition at line 184 of file mesh_smoother_vsmoother.h.

Referenced by smooth().

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

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

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

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

bool libMesh::VariationalMeshSmoother::_generate_data
private

Definition at line 186 of file mesh_smoother_vsmoother.h.

Referenced by set_generate_data(), and smooth().

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

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

const unsigned libMesh::VariationalMeshSmoother::_maxiter
private

Definition at line 181 of file mesh_smoother_vsmoother.h.

Referenced by smooth().

MetricType libMesh::VariationalMeshSmoother::_metric
private

Definition at line 183 of file mesh_smoother_vsmoother.h.

Referenced by set_metric(), and smooth().

const unsigned libMesh::VariationalMeshSmoother::_miniter
private

Definition at line 180 of file mesh_smoother_vsmoother.h.

Referenced by smooth().

const unsigned libMesh::VariationalMeshSmoother::_miniterBC
private

Definition at line 182 of file mesh_smoother_vsmoother.h.

Referenced by smooth().

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

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

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

const double libMesh::VariationalMeshSmoother::_percent_to_move
private

Dampening factor

Definition at line 159 of file mesh_smoother_vsmoother.h.

Referenced by writegr().

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: