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

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

#include <mesh_smoother_vsmoother.h>

Inheritance diagram for libMesh::VariationalMeshSmoother:
[legend]

Classes

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

Public Types

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

Public Member Functions

 VariationalMeshSmoother (UnstructuredMesh &mesh, double theta=0.5, unsigned miniter=2, unsigned maxiter=5, unsigned miniterBC=5)
 Simple constructor to use for smoothing purposes. More...
 
 VariationalMeshSmoother (UnstructuredMesh &mesh, std::vector< float > *adapt_data, 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. More...
 
 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. More...
 
virtual ~VariationalMeshSmoother ()
 Destructor. More...
 
virtual void smooth () libmesh_override
 Redefinition of the smooth function from the base class. More...
 
double smooth (unsigned int n_iterations)
 The actual smoothing function, gets called whenever the user specifies an actual number of smoothing iterations. More...
 
double distance_moved () const
 
void set_generate_data (bool b)
 Allow user to control whether the metric is generated from the initial mesh. More...
 
void set_metric (MetricType t)
 Allow user to control the smoothing metric used. More...
 

Protected Attributes

UnstructuredMesh_mesh
 

Private Member Functions

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

Detailed Description

This is an implementation of Larisa Branets' smoothing algorithms.

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

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

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

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

Author
Derek R. Gaston
Date
2006

Definition at line 63 of file mesh_smoother_vsmoother.h.

Member Enumeration Documentation

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

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

References _adapt_data, and std::min().

Referenced by adjust_adapt_data().

517 {
518  float min = 1.e30;
519 
520  for (std::size_t i=0; i<_adapt_data->size(); i++)
521  {
522  // Only positive (or zero) values in the error vector
523  libmesh_assert_greater_equal ((*_adapt_data)[i], 0.);
524  min = std::min (min, (*_adapt_data)[i]);
525  }
526 
527  // ErrorVectors are for positive values
528  libmesh_assert_greater_equal (min, 0.);
529 
530  return min;
531 }
std::vector< float > * _adapt_data
Vector for holding adaptive data.
long double min(long double a, double b)
void libMesh::VariationalMeshSmoother::adjust_adapt_data ( )
private

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

536 {
537  // For convenience
538  const UnstructuredMesh & aoe_mesh = *_area_of_interest;
539  std::vector<float> & adapt_data = *_adapt_data;
540 
541  float min = adapt_minimum();
542 
543  MeshBase::const_element_iterator el = _mesh.elements_begin();
544  const MeshBase::const_element_iterator end_el = _mesh.elements_end();
545 
546  MeshBase::const_element_iterator aoe_el = aoe_mesh.elements_begin();
547  const MeshBase::const_element_iterator aoe_end_el = aoe_mesh.elements_end();
548 
549  // Counter to keep track of which spot in adapt_data we're looking at
550  for (int i=0; el!=end_el; el++)
551  {
552  // Only do this for active elements
553  if (adapt_data[i])
554  {
555  Point centroid = (*el)->centroid();
556  bool in_aoe = false;
557 
558  // See if the elements centroid lies in the aoe mesh
559  // for (aoe_el=aoe_mesh.elements_begin(); aoe_el != aoe_end_el; ++aoe_el)
560  // {
561  if ((*aoe_el)->contains_point(centroid))
562  in_aoe = true;
563  // }
564 
565  // If the element is not in the area of interest... then set
566  // it's adaptivity value to the minimum.
567  if (!in_aoe)
568  adapt_data[i] = min;
569  }
570  i++;
571  }
572 }
UnstructuredMesh & _mesh
Definition: mesh_smoother.h:61
virtual element_iterator elements_begin()=0
Iterate over all the elements in the Mesh.
virtual element_iterator elements_end()=0
std::vector< float > * _adapt_data
Vector for holding adaptive data.
const UnstructuredMesh * _area_of_interest
Area of Interest Mesh.
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 833 of file mesh_smoother_vsmoother.C.

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

Referenced by full_smooth().

837 {
838  // evaluates given adaptive function on the provided mesh
839  if (adp < 0)
840  {
841  for (dof_id_type i=0; i<_n_cells; i++)
842  {
843  double
844  x = 0.,
845  y = 0.,
846  z = 0.;
847  int nvert = 0;
848 
849  while (cells[i][nvert] >= 0)
850  {
851  x += R[cells[i][nvert]][0];
852  y += R[cells[i][nvert]][1];
853  if (_dim == 3)
854  z += R[cells[i][nvert]][2];
855  nvert++;
856  }
857 
858  // adaptive function, cell based
859  afun[i] = 5.*(x+y+z);
860  }
861  }
862 
863  else if (adp > 0)
864  {
865  for (dof_id_type i=0; i<_n_nodes; i++)
866  {
867  double z = 0;
868  for (unsigned j=0; j<_dim; j++)
869  z += R[i][j];
870 
871  // adaptive function, node based
872  afun[i] = 5*std::sin(R[i][0]);
873  }
874  }
875 }
dof_id_type _n_cells
The number of active elements in the Mesh at the time of smoothing.
dof_id_type _n_nodes
The number of nodes in the Mesh at the time of smoothing.
PetscErrorCode Vec x
uint8_t dof_id_type
Definition: id_types.h:64
const unsigned _dim
Smoother control variables.
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 3301 of file mesh_smoother_vsmoother.C.

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

Referenced by localP().

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

Definition at line 625 of file mesh_smoother_vsmoother.C.

References _dim.

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

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

Definition at line 137 of file mesh_smoother_vsmoother.h.

References _distance.

137 { return _distance; }
double _distance
Max distance of the last set of movement.
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 880 of file mesh_smoother_vsmoother.C.

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

Referenced by smooth().

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

References libMesh::x.

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

Definition at line 614 of file mesh_smoother_vsmoother.C.

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

618 {
619  return x1*y2 - x2*y1;
620 }
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 599 of file mesh_smoother_vsmoother.C.

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

608 {
609  return x1*(y2*z3 - y3*z2) + y1*(z2*x3 - z3*x2) + z1*(x2*y3 - x3*y2);
610 }
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 2970 of file mesh_smoother_vsmoother.C.

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

Referenced by minJ(), and minJ_BC().

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

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

Referenced by full_smooth().

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

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

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

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

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

Referenced by full_smooth().

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

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

Referenced by full_smooth().

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

References _logfile, and libMesh::x.

Referenced by solver().

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

References _logfile.

Referenced by solver().

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

Definition at line 577 of file mesh_smoother_vsmoother.C.

References _adapt_data, _area_of_interest, and adjust_adapt_data().

Referenced by full_smooth().

578 {
579  std::vector<float> & adapt_data = *_adapt_data;
580 
581  if (_area_of_interest)
583 
584  std::size_t m = adapt_data.size();
585 
586  std::size_t j =0 ;
587  for (std::size_t i=0; i<m; i++)
588  if (adapt_data[i] != 0)
589  {
590  afun[j] = adapt_data[i];
591  j++;
592  }
593 
594  return 0;
595 }
std::vector< float > * _adapt_data
Vector for holding adaptive data.
const UnstructuredMesh * _area_of_interest
Area of Interest Mesh.
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 268 of file mesh_smoother_vsmoother.C.

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

Referenced by metr_data_gen(), and smooth().

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

Definition at line 494 of file mesh_smoother_vsmoother.C.

References _dim, and _n_cells.

Referenced by smooth().

496 {
497  std::ifstream infile(name.c_str());
498  std::string dummy;
499 
500  for (dof_id_type i=0; i<_n_cells; i++)
501  for (unsigned j=0; j<_dim; j++)
502  {
503  for (unsigned k=0; k<_dim; k++)
504  infile >> H[i][j][k];
505 
506  // Read to end of line and discard
507  std::getline(infile, dummy);
508  }
509 
510  return 0;
511 }
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:39
dof_id_type _n_cells
The number of active elements in the Mesh at the time of smoothing.
uint8_t dof_id_type
Definition: id_types.h:64
const unsigned _dim
Smoother control variables.
void libMesh::VariationalMeshSmoother::set_generate_data ( bool  b)

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)

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

Redefinition of the smooth function from the base class.

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

Implements libMesh::MeshSmoother.

Definition at line 125 of file mesh_smoother_vsmoother.h.

References _distance, and smooth().

Referenced by smooth().

125 { _distance = this->smooth(1); }
double _distance
Max distance of the last set of movement.
virtual void smooth() libmesh_override
Redefinition of the smooth function from the base class.
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
Map for 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
dof_id_type _n_hanging_edges
The number of hanging node edges in the Mesh at the time of smoothing.
std::ofstream _logfile
All output (including debugging) is sent to the _logfile.
UnstructuredMesh & _mesh
Definition: mesh_smoother.h:61
dof_id_type _n_cells
The number of active elements in the Mesh at the time of smoothing.
dof_id_type _n_nodes
The number of nodes in the Mesh at the time of smoothing.
int readgr(Array2D< 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)
Given a mesh hanging_nodes will be filled with an associative array keyed off the global id of all th...
Definition: mesh_tools.C:873
virtual dof_id_type n_nodes() const =0
int writegr(const Array2D< double > &R)
const unsigned _dim
Smoother control variables.
double _dist_norm
Records a relative "distance moved".
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 3611 of file mesh_smoother_vsmoother.C.

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

Referenced by minJ().

3620 {
3621  _logfile << "Beginning Solve " << n << std::endl;
3622 
3623  // Check parameters
3624  int info = pcg_par_check(n, ia, ja, a, eps, maxite, msglev);
3625  if (info != 0)
3626  return info;
3627 
3628  // PJ preconditioner construction
3629  std::vector<double> u(n);
3630  for (int i=0; i<n; i++)
3631  u[i] = 1./a[ia[i]];
3632 
3633  // PCG iterations
3634  std::vector<double> r(n), p(n), z(n);
3635  info = pcg_ic0(n, ia, ja, a, u, x, b, r, p, z, eps, maxite, msglev);
3636 
3637  // Mat sparse_global;
3638  // MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,n,n,ia,ja,a,&sparse_global);
3639 
3640  _logfile << "Finished Solve " << n << std::endl;
3641 
3642  return info;
3643 }
std::ofstream _logfile
All output (including debugging) is sent to the _logfile.
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 3393 of file mesh_smoother_vsmoother.C.

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

Referenced by localP().

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

Referenced by smooth().

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

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

UnstructuredMesh& libMesh::MeshSmoother::_mesh
protectedinherited
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: