libMesh::VariationalMeshSmoother Class Reference

#include <mesh_smoother_vsmoother.h>

Inheritance diagram for libMesh::VariationalMeshSmoother:

Public Types
enum	metric_type { uniform = 1, volumetric = 2, directional = 3 }
enum	adapt_type { cell = -1, none = 0, node = 1 }
Public Member Functions
	VariationalMeshSmoother (UnstructuredMesh &mesh, const double theta=0.5, const uint miniter=2, const uint maxiter=5, const uint miniterBC=5)
	VariationalMeshSmoother (UnstructuredMesh &mesh, std::vector< float > *adapt_data, const double theta=0.5, const uint miniter=2, const uint maxiter=5, const uint miniterBC=5, const double percent_to_move=1)
	VariationalMeshSmoother (UnstructuredMesh &mesh, const UnstructuredMesh area_of_interest, std::vector< float > adapt_data, const double theta=0.5, const uint miniter=2, const uint maxiter=5, const uint miniterBC=5, const double percent_to_move=1)
virtual	~VariationalMeshSmoother ()
virtual void	smooth ()
double	smooth (unsigned int n_iterations)
double	distanceMoved () const
	Member function `distanceMoved`.
Protected Attributes
UnstructuredMesh &	_mesh
Private Member Functions
void	adjust_adapt_data ()
float	adapt_minimum () const
LPDOUBLE	alloc_d_n1 (int m1)
LPINT	alloc_i_n1 (int m1)
LPLPINT	alloc_i_n1_n2 (int m1, int)
LPLPDOUBLE	alloc_d_n1_n2 (int m1, int)
LPLPLPDOUBLE	alloc_d_n1_n2_n3 (int m1, int, int)
int	writegr (int n, int N, LPLPDOUBLE R, LPINT mask, int ncells, LPLPINT cells, LPINT mcells, int nedges, LPINT edges, LPINT hnodes, const char grid[], int me, const char grid_old[], FILE *sout)
int	readgr (int n, int N, LPLPDOUBLE R, LPINT mask, int ncells, LPLPINT cells, LPINT mcells, int nedges, LPINT edges, LPINT hnodes, FILE *sout)
int	readmetr (char name, LPLPLPDOUBLE H, int ncells, int n, FILE sout)
int	read_adp (LPDOUBLE afun, char adap, FILE sout)
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 (int n, LPLPDOUBLE Q, int nvert, LPDOUBLE K, LPLPDOUBLE H, int me)
void	adp_renew (int n, int N, LPLPDOUBLE R, int ncells, LPLPINT cells, LPDOUBLE afun, int adp, FILE *sout)
void	full_smooth (int n, int N, LPLPDOUBLE R, LPINT mask, int ncells, LPLPINT cells, LPINT mcells, int nedges, LPINT edges, LPINT hnodes, double w, LPINT iter, int me, LPLPLPDOUBLE H, int adp, char adap, int gr, FILE sout)
double	maxE (int n, int N, LPLPDOUBLE R, int ncells, LPLPINT cells, LPINT mcells, int me, LPLPLPDOUBLE H, double v, double epsilon, double w, LPDOUBLE Gamma, double qmin, FILE sout)
double	minq (int n, int N, LPLPDOUBLE R, int ncells, LPLPINT cells, LPINT mcells, int me, LPLPLPDOUBLE H, double vol, double Vmin, FILE *sout)
double	minJ (int n, int N, LPLPDOUBLE R, LPINT mask, int ncells, LPLPINT cells, LPINT mcells, double epsilon, double w, int me, LPLPLPDOUBLE H, double vol, int nedges, LPINT edges, LPINT hnodes, int msglev, double Vmin, double emax, double qmin, int adp, LPDOUBLE afun, FILE sout)
double	minJ_BC (int N, LPLPDOUBLE R, LPINT mask, int ncells, LPLPINT cells, LPINT mcells, double epsilon, double w, int me, LPLPLPDOUBLE H, double vol, int msglev, double Vmin, double emax, double qmin, int adp, LPDOUBLE afun, int NCN, FILE sout)
double	localP (int n, LPLPLPDOUBLE W, LPLPDOUBLE F, LPLPDOUBLE R, LPINT cell, LPINT mask, double epsilon, double w, int nvert, LPLPDOUBLE H, int me, double vol, int f, double Vmin, double qmin, int adp, LPDOUBLE afun, LPDOUBLE Gloc, FILE *sout)
double	avertex (int n, LPDOUBLE afun, LPDOUBLE G, LPLPDOUBLE R, LPINT cell, int nvert, int adp, FILE *sout)
double	vertex (int n, LPLPLPDOUBLE W, LPLPDOUBLE F, LPLPDOUBLE R, LPINT cell, double epsilon, double w, int nvert, LPDOUBLE K, LPLPDOUBLE H, int me, double vol, int f, double Vmin, int adp, LPDOUBLE G, double sigma, FILE sout)
void	metr_data_gen (char grid[], char metr[], int n, int me, FILE *sout)
int	solver (int n, LPINT ia, LPINT ja, LPDOUBLE a, LPDOUBLE x, LPDOUBLE b, double eps, int maxite, int msglev, FILE *sout)
int	pcg_ic0 (int n, LPINT ia, LPINT ja, LPDOUBLE a, LPDOUBLE u, LPDOUBLE x, LPDOUBLE b, LPDOUBLE r, LPDOUBLE p, LPDOUBLE z, double eps, int maxite, int msglev, FILE *sout)
int	pcg_par_check (int n, LPINT ia, LPINT ja, LPDOUBLE a, double eps, int maxite, int msglev, FILE *sout)
void	gener (char grid[], int n, FILE *sout)
void	local_sweep (int n, int N, LPLPDOUBLE R, LPINT mask, int ncells, LPLPINT cells, LPINT mcells, int nedges, LPINT edges, LPINT hnodes, double w, LPINT iter, int me, LPLPLPDOUBLE H, int adp, int OPT, FILE *sout)
double	minJ_l (int n, int N, LPLPDOUBLE R, LPINT mask, int ncells, LPLPINT cells, LPINT mcells, double epsilon, double w, int me, LPLPLPDOUBLE H, double vol, int nedges, LPINT edges, LPINT hnodes, int msglev, double Vmin, double emax, double qmin, int adp, LPDOUBLE afun, FILE sout)
Private Attributes
double	_distance
const double	_percent_to_move
double	_dist_norm
std::map< dof_id_type, std::vector< dof_id_type > >	_hanging_nodes
std::vector< float > *	_adapt_data
const uint	_dim
const uint	_miniter
const uint	_maxiter
const uint	_miniterBC
const metric_type	_metric
const adapt_type	_adaptive_func
const double	_theta
const bool	_generate_data
const UnstructuredMesh *	_area_of_interest

Detailed Description

This is an implementation of Larisa Branets smoothing algorithms. The initial implementation was done by her, the adaptation to libmesh was completed by Derek Gaston.

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.V. 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 71 of file mesh_smoother_vsmoother.h.

Member Enumeration Documentation

enum libMesh::VariationalMeshSmoother::adapt_type

Enumerator:

cell
none
node

Definition at line 148 of file mesh_smoother_vsmoother.h.

00149   {
00150     cell=-1,
00151     none=0,
00152     node=1
00153   };

enum libMesh::VariationalMeshSmoother::metric_type

Enumerator:

uniform
volumetric
directional

Definition at line 141 of file mesh_smoother_vsmoother.h.

00142   {
00143     uniform=1,
00144     volumetric=2,
00145     directional=3
00146   };

Constructor & Destructor Documentation

libMesh::VariationalMeshSmoother::VariationalMeshSmoother	(	UnstructuredMesh &	mesh,
		const double	theta = `0.5`,
		const uint	miniter = `2`,
		const uint	maxiter = `5`,
		const uint	miniterBC = `5`
	)			`[inline]`

Simple constructor to use for smoothing purposes

Definition at line 78 of file mesh_smoother_vsmoother.h.

00080     :MeshSmoother(mesh),
00081      _percent_to_move(1),
00082      _adapt_data(NULL),
00083      _dim(mesh.mesh_dimension()),
00084      _miniter(miniter),
00085      _maxiter(maxiter),
00086      _miniterBC(miniterBC),
00087 
00088      _metric(uniform),
00089      _adaptive_func(none),
00090      _theta(theta),
00091      _generate_data(false),
00092 
00093      _area_of_interest(NULL)
00094   {}

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

Slightly more complicated constructor for mesh redistribution based on adapt_data

Definition at line 99 of file mesh_smoother_vsmoother.h.

00102     :MeshSmoother(mesh),
00103      _percent_to_move(percent_to_move),
00104      _adapt_data(adapt_data),
00105      _dim(mesh.mesh_dimension()),
00106      _miniter(miniter),
00107      _maxiter(maxiter),
00108      _miniterBC(miniterBC),
00109 
00110      _metric(uniform),
00111      _adaptive_func(cell),
00112      _theta(theta),
00113      _generate_data(false),
00114 
00115      _area_of_interest(NULL)
00116   {}

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

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

Definition at line 122 of file mesh_smoother_vsmoother.h.

00125     :MeshSmoother(mesh),
00126      _percent_to_move(percent_to_move),
00127      _adapt_data(adapt_data),
00128      _dim(mesh.mesh_dimension()),
00129      _miniter(miniter),
00130      _maxiter(maxiter),
00131      _miniterBC(miniterBC),
00132 
00133      _metric(uniform),
00134      _adaptive_func(cell),
00135      _theta(theta),
00136      _generate_data(false),
00137 
00138      _area_of_interest(area_of_interest)
00139   {}

virtual libMesh::VariationalMeshSmoother::~VariationalMeshSmoother ( ) [inline, virtual]

Destructor.

Definition at line 158 of file mesh_smoother_vsmoother.h.

00158 {}

Member Function Documentation

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

Definition at line 520 of file mesh_smoother_vsmoother.C.

References _adapt_data, and std::min().

Referenced by adjust_adapt_data().

00521 {
00522   std::vector<float>& adapt_data = *_adapt_data;
00523 
00524   const std::size_t n = adapt_data.size();
00525   float min = 1.e30;
00526 
00527   for (unsigned int i=0; i<n; i++)
00528   {
00529       // Only positive (or zero) values in the error vector
00530     libmesh_assert_greater_equal (adapt_data[i], 0.);
00531     min = std::min (min, adapt_data[i]);
00532   }
00533 
00534   // ErrorVectors are for positive values
00535   libmesh_assert_greater_equal (min, 0.);
00536 
00537   return min;
00538 }

void libMesh::VariationalMeshSmoother::adjust_adapt_data ( ) [private]

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

00541 {
00542   //For convenience
00543   const UnstructuredMesh & aoe_mesh=*_area_of_interest;
00544   std::vector<float>& adapt_data = *_adapt_data;
00545 
00546   float min=adapt_minimum();
00547 
00548   MeshBase::const_element_iterator       el     = _mesh.elements_begin();
00549   const MeshBase::const_element_iterator end_el = _mesh.elements_end();
00550 
00551   MeshBase::const_element_iterator       aoe_el     = aoe_mesh.elements_begin();
00552   const MeshBase::const_element_iterator aoe_end_el = aoe_mesh.elements_end();
00553 
00554   //Counter to keep track of which spot in adapt_data we're looking at
00555   int i=0;
00556 
00557   for(; el != end_el; el++)
00558   {
00559     //Only do this for active elements
00560     if(adapt_data[i])
00561     {
00562       Point centroid=(*el)->centroid();
00563       bool in_aoe=false;
00564 
00565       //See if the elements centroid lies in the aoe mesh
00566       //for(aoe_el=aoe_mesh.elements_begin(); aoe_el != aoe_end_el; ++aoe_el)
00567       //{
00568         if((*aoe_el)->contains_point(centroid))
00569           in_aoe=true;
00570       //}
00571 
00572       //If the element is not in the area of interest... then set
00573       //it's adaptivity value to the minimum.
00574       if(!in_aoe)
00575         adapt_data[i]=min;
00576     }
00577     i++;
00578   }
00579 }

void libMesh::VariationalMeshSmoother::adp_renew	(	int	n,
		int	N,
		LPLPDOUBLE	R,
		int	ncells,
		LPLPINT	cells,
		LPDOUBLE	afun,
		int	adp,
		FILE *	sout
	)			`[private]`

Specify adaptive function

Definition at line 761 of file mesh_smoother_vsmoother.C.

00762 {
00763   int i, j, nvert;
00764   double x, y, z;
00765 
00766   //evaluates given adaptive function on the provided mesh
00767   if(adp<0){
00768     for(i=0;i<ncells;i++)
00769     {
00770       x=y=z=0;
00771       nvert=0;
00772       while(cells[i][nvert]>=0)
00773       {
00774         x+=R[cells[i][nvert]][0];
00775         y+=R[cells[i][nvert]][1];
00776         if(n==3)
00777           z+=R[cells[i][nvert]][2];
00778         nvert++;
00779       }
00780       afun[i]=5*(x+y+z); //adaptive function, cell based
00781     }
00782   }
00783   if(adp>0)
00784   {
00785     for(i=0;i<N;i++)
00786     {
00787       z=0; for(j=0;j<n;j++) z+=R[i][j];
00788       afun[i]=5*sin(R[i][0]); //adaptive function, node based
00789     }
00790   }
00791 
00792   return;
00793 }

LPDOUBLE libMesh::VariationalMeshSmoother::alloc_d_n1 ( int m1 ) [inline, private]

Imported stuff: ..........

Definition at line 231 of file mesh_smoother_vsmoother.h.

Referenced by avertex(), basisA(), maxE(), metr_data_gen(), minJ(), minJ_BC(), minq(), smooth(), solver(), and vertex().

00231 { return((LPDOUBLE)malloc(m1*sizeof(double))); }

LPLPDOUBLE libMesh::VariationalMeshSmoother::alloc_d_n1_n2	(	int	m1,
		int
	)			`[inline, private]`

Definition at line 234 of file mesh_smoother_vsmoother.h.

Referenced by avertex(), basisA(), maxE(), metr_data_gen(), minJ(), minJ_BC(), minq(), smooth(), and vertex().

00234 { return((LPLPDOUBLE)malloc(m1*sizeof(LPDOUBLE))); }

LPLPLPDOUBLE libMesh::VariationalMeshSmoother::alloc_d_n1_n2_n3	(	int	m1,
		int	,
		int
	)			`[inline, private]`

Definition at line 235 of file mesh_smoother_vsmoother.h.

Referenced by minJ(), minJ_BC(), smooth(), and vertex().

00235 { return((LPLPLPDOUBLE)malloc(m1*sizeof(LPLPDOUBLE))); }

LPINT libMesh::VariationalMeshSmoother::alloc_i_n1 ( int m1 ) [inline, private]

Definition at line 232 of file mesh_smoother_vsmoother.h.

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

00232 { return((LPINT)malloc(m1*sizeof(int))); }

LPLPINT libMesh::VariationalMeshSmoother::alloc_i_n1_n2	(	int	m1,
		int
	)			`[inline, private]`

Definition at line 233 of file mesh_smoother_vsmoother.h.

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

00233 { return((LPLPINT)malloc(m1*sizeof(LPINT))); }

double libMesh::VariationalMeshSmoother::avertex	(	int	n,
		LPDOUBLE	afun,
		LPDOUBLE	G,
		LPLPDOUBLE	R,
		LPINT	cell_in,
		int	nvert,
		int	adp,
		FILE *	sout
	)			`[private]`

avertex - assembly of adaptivity metric on a cell

Definition at line 2112 of file mesh_smoother_vsmoother.C.

References alloc_d_n1(), alloc_d_n1_n2(), basisA(), jac2(), and jac3().

Referenced by localP().

02113 {
02114   LPLPDOUBLE Q;
02115   LPDOUBLE K;
02116   double a1[3], a2[3], a3[3], qu[3];
02117   int i,j;
02118   double det, g, df0, df1, df2;
02119   K=alloc_d_n1(8);
02120   Q = alloc_d_n1_n2(3, nvert);
02121   for(i=0;i<n;i++) Q[i]=alloc_d_n1(nvert);
02122 
02123   for(i=0;i<8;i++) K[i]=0.5; //cell center
02124 
02125   basisA(n,Q,nvert,K,Q,1);
02126   for(i=0;i<n;i++){a1[i]=0; a2[i]=0; a3[i]=0; qu[i]=0;
02127      for(j=0;j<nvert;j++){
02128         a1[i]+=Q[i][j]*R[cell_in[j]][0];
02129         a2[i]+=Q[i][j]*R[cell_in[j]][1];
02130                 if(n==3) a3[i]+=Q[i][j]*R[cell_in[j]][2];
02131                 qu[i]+=Q[i][j]*afun[cell_in[j]];
02132      }
02133   }
02134   if(n==2) det=jac2(a1[0],a1[1],a2[0],a2[1]); else det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
02135   if(det!=0){
02136           if(n==2){
02137           df0=jac2(qu[0],qu[1],a2[0],a2[1])/det;
02138           df1=jac2(a1[0],a1[1],qu[0],qu[1])/det;
02139                   g=(1+df0*df0+df1*df1);
02140                   if(adp==2) { //directional adaptation G=diag(g_i)
02141                           G[0]=1+df0*df0; G[1]=1+df1*df1;
02142                   } else {
02143                           for(i=0;i<n;i++) G[i]=g; //simple adaptation G=gI
02144                   }
02145           } else {
02146                   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;
02147           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;
02148           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;
02149                   g=(1+df0*df0+df1*df1+df2*df2);
02150           if(adp==2){//directional adaptation G=diag(g_i)
02151                           G[0]=1+df0*df0; G[1]=1+df1*df1; G[2]=1+df2*df2;
02152                   } else {
02153                           for(i=0;i<n;i++) G[i]=g; //simple adaptation G=gI
02154                   }
02155           }
02156   } else {g=1.0; for(i=0;i<n;i++) G[i]=g;}
02157 
02158 
02159  for(i=0;i<n;i++) free(Q[i]);
02160  free(Q);
02161 
02162  return(g);
02163 
02164 }

int libMesh::VariationalMeshSmoother::basisA	(	int	n,
		LPLPDOUBLE	Q,
		int	nvert,
		LPDOUBLE	K,
		LPLPDOUBLE	H,
		int	me
	)			`[private]`

BasisA determines matrix H^(-T)Q on one Jacobian matrix

Definition at line 623 of file mesh_smoother_vsmoother.C.

References alloc_d_n1(), and alloc_d_n1_n2().

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

00625 00626 00627 00628 00629 00630 00631 00632 00633 00634 00635 00636 00637 00638 00639 00640     } 00641 00642 00643 00644 00645 00646 00647 00648 00649 00650     } 00651 00652 00653 00654 00655 00656 00657 00658 00659 00660 00661 00662 00663 00664     } 00665  } 00666 00667 00668 00669 00670 00671 00672 00673 00674 00675 00676 00677 00678 00679 00680 00681 00682 00683 00684 00685 00686 00687 00688 00689 00690 00691 00692      } 00693 00694 00695 00696 00697 00698 00699 00700 00701 00702 00703 00704      } 00705 00706 00707 00708 00709 00710 00711 00712 00713      } 00714 00715 00716 00717 00718 00719 00720 00721 00722 00723 00724 00725 00726 00727 00728 00729 00730 00731 00732 00733 00734 00735 00736 00737      } 00738  } 00739 00740 00741 00742 00743 00744 00745 00746 00747 00748        } 00749     } 00750  } 00751 00752 00753  free(U); 00754 00755 }

class="fragment">00624 { LPLPDOUBLE U; int i,j,k; U=alloc_d_n1_n2(n,nvert); for(i=0;i<n;i++) U[i]=alloc_d_n1(nvert); if(n==2){ if(nvert==4){//quad U[0][0]=(-(1-K[1])); U[0][1]=( (1-K[1])); U[0][2]=(- K[1]); U[0][3]=( K[1]); U[1][0]=(-(1-K[0])); U[1][1]=(- K[0]); U[1][2]=( (1-K[0])); U[1][3]=( K[0]); if(nvert==3){//tri /*---for right target triangle--- U[0][0]=-1.0; U[1][0]=-1.0; U[0][1]=1.0; U[1][1]=0.0; U[0][2]=0.0; U[1][2]=1.0; -----for regular triangle---*/ U[0][0]=-1.0; U[1][0]=-1.0/sqrt(3.0); U[0][1]= 1.0; U[1][1]=-1.0/sqrt(3.0); U[0][2]= 0; U[1][2]= 2.0/sqrt(3.0); if(nvert==6){/*--curved triangle*/ U[0][0]=-3*(1-K[0]-K[1]*0.5)+(K[0]-0.5*K[1])+K[1]; U[0][1]=-(1-K[0]-K[1]*0.5)+(K[0]-0.5*K[1])*3-K[1]; U[0][2]=0; U[0][3]=K[1]*4; U[0][4]=-4*K[1]; U[0][5]=4*(1-K[0]-K[1]*0.5)-4*(K[0]-0.5*K[1]); U[1][0]=(1-K[2]-K[3]*0.5)*(-1.5)+(K[2]-0.5*K[3])*(0.5)+K[3]*(0.5); U[1][1]=(1-K[2]-K[3]*0.5)*(0.5)+(K[2]-0.5*K[3])*(-1.5)+K[3]*(0.5); U[1][2]=(1-K[2]-K[3]*0.5)*(-1)+(K[2]-0.5*K[3])*(-1)+K[3]*(3); U[1][3]=(K[2]-0.5*K[3])*(4.0)+K[3]*(-2.0); U[1][4]=(1-K[2]-K[3]*0.5)*(4.0)+K[3]*(-2.0); U[1][5]=(1-K[2]-K[3]*0.5)*(-2.0)+(K[2]-0.5*K[3])*(-2.0); if(n==3){ if(nvert==8){//hex U[0][0]=(-(1-K[0])*(1-K[1])); U[0][1]=( (1-K[0])*(1-K[1])); U[0][2]=(- K[0]*(1-K[1])); U[0][3]=( K[0]*(1-K[1])); U[0][4]=(-(1-K[0])* K[1]); U[0][5]=( (1-K[0])* K[1]); U[0][6]=(- K[0]* K[1]); U[0][7]=( K[0]* K[1]); U[1][0]=(-(1-K[2])*(1-K[3])); U[1][1]=(- K[2]*(1-K[3])); U[1][2]=( (1-K[2])*(1-K[3])); U[1][3]=( K[2]*(1-K[3])); U[1][4]=(-(1-K[2])* K[3]); U[1][5]=(- K[2]* K[3]); U[1][6]=( (1-K[2])* K[3]); U[1][7]=( K[2]* K[3]); U[2][0]=(-(1-K[4])*(1-K[5])); U[2][1]=(- K[4]*(1-K[5])); U[2][2]=(-(1-K[4])* K[5]); U[2][3]=(- K[4]* K[5]); U[2][4]=( (1-K[4])*(1-K[5])); U[2][5]=( K[4]*(1-K[5])); U[2][6]=( (1-K[4])* K[5]); U[2][7]=( K[4]* K[5]); if(nvert==4){//linear tetr /*---for regular reference tetrahedron--------*/ U[0][0]=-1; U[1][0]=-1.0/sqrt(3.0); U[2][0]=-1.0/sqrt(6.0); U[0][1]=1; U[1][1]=-1.0/sqrt(3.0); U[2][1]=-1.0/sqrt(6.0); U[0][2]=0; U[1][2]=2.0/sqrt(3.0); U[2][2]=-1.0/sqrt(6.0); U[0][3]=0; U[1][3]=0; U[2][3]=3.0/sqrt(6.0); /*---for right corner reference tetrahedron---- U[0][0]=-1; U[1][0]=-1; U[2][0]=-1; U[0][1]=1; U[1][1]=0; U[2][1]=0; U[0][2]=0; U[1][2]=1; U[2][2]=0; U[0][3]=0; U[1][3]=0; U[2][3]=1;*/ if(nvert==6){//prism /* for regular triangle in the prism base */ U[0][0]=-1+K[0]; U[1][0]=(-1+K[1])/2.0; U[2][0]=-1+K[2]+K[3]/2.0; U[0][1]=1-K[0]; U[1][1]=(-1+K[1])/2.0; U[2][1]=-K[2]+K[3]/2.0; U[0][2]=0; U[1][2]=(1-K[1]); U[2][2]=-K[3]; U[0][3]=-K[0]; U[1][3]=(-K[1])/2.0; U[2][3]=1-K[2]-K[3]/2.0; U[0][4]=K[0]; U[1][4]=(-K[1])/2.0; U[2][4]=K[2]-K[3]/2.0; U[0][5]=0; U[1][5]=K[1]; U[2][5]=K[3]; if(nvert==10){//quad tetr U[0][0]=(1-K[0]-K[1]/2-K[2]/3)*(-3)+(K[0]-K[1]/2-K[2]/3)+(K[1]-K[2]/3)+K[2]; U[0][1]=(1-K[0]-K[1]/2-K[2]/3)*(-1)+(K[0]-K[1]/2-K[2]/3)*3-(K[1]-K[2]/3)-K[2]; U[0][2]=0; U[0][3]=0; U[0][4]=4*(K[1]-K[2]/3); U[0][5]=-4*(K[1]-K[2]/3); U[0][6]=4*(1-K[0]-K[1]/2-K[2]/3)-4*(K[0]-K[1]/2-K[2]/3); U[0][7]=4*K[2]; U[0][8]=0; U[0][9]=-4*K[2]; U[1][0]=(1-K[3]-K[4]/2-K[5]/3)*(-3/2)+(K[3]-K[4]/2-K[5]/3)/2+(K[4]-K[5]/3)/2+K[5]/2; U[1][1]=(1-K[3]-K[4]/2-K[5]/3)/2+(K[3]-K[4]/2-K[5]/3)*(-3/2)+(K[4]-K[5]/3)/2+K[5]/2; U[1][2]=-(1-K[3]-K[4]/2-K[5]/3)-(K[3]-K[4]/2-K[5]/3)+3*(K[4]-K[5]/3)-K[5]; U[1][3]=0; U[1][4]=4*(K[3]-K[4]/2-K[5]/3)-2*(K[4]-K[5]/3); U[1][5]=4*(1-K[3]-K[4]/2-K[5]/3)-2*(K[4]-K[5]/3); U[1][6]=-2*(1-K[3]-K[4]/2-K[5]/3)-2*(K[3]-K[4]/2-K[5]/3); U[1][7]=-2*K[5]; U[1][8]=4*K[5]; U[1][9]=-2*K[5]; U[2][0]=-(1-K[6]-K[7]/2-K[8]/3)+(K[6]-K[7]/2-K[8]/3)/3+(K[7]-K[8]/3)/3+K[8]/3; U[2][1]=(1-K[6]-K[7]/2-K[8]/3)/3-(K[6]-K[7]/2-K[8]/3)+(K[7]-K[8]/3)/3+K[8]/3; U[2][2]=(1-K[6]-K[7]/2-K[8]/3)/3+(K[6]-K[7]/2-K[8]/3)/3-(K[7]-K[8]/3)+K[8]/3; U[2][3]=-(1-K[6]-K[7]/2-K[8]/3)-(K[6]-K[7]/2-K[8]/3)-(K[7]-K[8]/3)+3*K[8]; U[2][4]=-4*(K[6]-K[7]/2-K[8]/3)/3-4*(K[7]-K[8]/3)/3; U[2][5]=-4*(1-K[6]-K[7]/2-K[8]/3)/3-4*(K[7]-K[8]/3)/3; U[2][6]=-4*(1-K[6]-K[7]/2-K[8]/3)/3-4*(K[6]-K[7]/2-K[8]/3)/3; U[2][7]=4*(K[6]-K[7]/2-K[8]/3)-4*K[8]/3; U[2][8]=4*(K[7]-K[8]/3)-4*K[8]/3; U[2][9]=4*(1-K[6]-K[7]/2-K[8]/3)-4*K[8]/3; if(me==1){ for(i=0;i<n;i++) for(j=0;j<nvert;j++) Q[i][j]=U[i][j]; } else { for(i=0;i<n;i++){ for(j=0;j<nvert;j++){ Q[i][j]=0; for(k=0;k<n;k++) Q[i][j]+=H[k][i]*U[k][j]; for(i=0;i<n;i++) free(U[i]); return 0;

double libMesh::VariationalMeshSmoother::distanceMoved ( ) const [inline]

Member function distanceMoved.

Returns:: a double max distance a node moved during the last smooth.

Definition at line 180 of file mesh_smoother_vsmoother.h.

References _distance.

00180 {return _distance;}

void libMesh::VariationalMeshSmoother::full_smooth	(	int	n,
		int	N,
		LPLPDOUBLE	R,
		LPINT	mask,
		int	ncells,
		LPLPINT	cells,
		LPINT	mcells,
		int	nedges,
		LPINT	edges,
		LPINT	hnodes,
		double	w,
		LPINT	iter,
		int	me,
		LPLPLPDOUBLE	H,
		int	adp,
		char *	adap,
		int	gr,
		FILE *	sout
	)			`[private]`

Referenced by smooth().

void libMesh::VariationalMeshSmoother::gener	(	char	grid[],
		int	n,
		FILE *	sout
	)			`[private]`

Sample mesh file generation

Definition at line 2586 of file mesh_smoother_vsmoother.C.

02587 {
02588         FILE *stream;
02589         int i, j, k, n1=3, N, ncells, mask, ii, jj, kk, nc;
02590         double x;
02591 
02592     N=1; ncells=1;
02593         for(i=0;i<n;i++){N*=n1; ncells*=(n1-1);}
02594         x=1.0/(double)(n1-1);
02595 
02596         stream=fopen(grid,"w+");
02597 
02598         fprintf(stream,"%d \n%d \n%d \n0 \n",n,N,ncells);
02599 
02600         for(i=0;i<N;i++){//node coordinates
02601                 k=i; mask=0;
02602                 for(j=0;j<n;j++){
02603                     ii=k%n1;
02604                     if((ii==0)||(ii==n1-1)) mask=1;
02605                     //if((i==N/2)&&(j==1)) fprintf(stream,"%e ",(double)ii*x+x/2.0); else
02606                     fprintf(stream,"%e ",(double)ii*x);
02607                     k/=n1;
02608                 }
02609         fprintf(stream,"%d \n",mask);
02610     }
02611     for(i=0;i<ncells;i++){//cell connectivity
02612             nc=i; ii=nc%(n1-1); nc/=(n1-1); jj=nc%(n1-1); kk=nc/(n1-1);
02613                 if(n==2) fprintf(stream,"%d %d %d %d ",ii+n1*jj,ii+1+jj*n1,ii+(jj+1)*n1,ii+1+(jj+1)*n1);
02614                 if(n==3) fprintf(stream,"%d %d %d %d %d %d %d %d ",ii+n1*jj+n1*n1*kk,ii+1+jj*n1+n1*n1*kk,ii+(jj+1)*n1+n1*n1*kk,ii+1+(jj+1)*n1+n1*n1*kk,
02615                                      ii+n1*jj+n1*n1*(kk+1),ii+1+jj*n1+n1*n1*(kk+1),ii+(jj+1)*n1+n1*n1*(kk+1),ii+1+(jj+1)*n1+n1*n1*(kk+1));
02616             fprintf(stream,"-1 -1 0 \n");
02617         }
02618         fclose(stream);
02619 
02620         return;
02621 }

double libMesh::VariationalMeshSmoother::jac2	(	double	x1,
		double	y1,
		double	x2,
		double	y2
	)			`[private]`

Definition at line 615 of file mesh_smoother_vsmoother.C.

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

00616 {
00617   return( x1*y2-x2*y1 );
00618 }

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

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

00611 {
00612   return( x1*(y2*z3 - y3*z2) + y1*(z2*x3 - z3*x2) + z1*(x2*y3 - x3*y2) );
00613 }

void libMesh::VariationalMeshSmoother::local_sweep	(	int	n,
		int	N,
		LPLPDOUBLE	R,
		LPINT	mask,
		int	ncells,
		LPLPINT	cells,
		LPINT	mcells,
		int	nedges,
		LPINT	edges,
		LPINT	hnodes,
		double	w,
		LPINT	iter,
		int	me,
		LPLPLPDOUBLE	H,
		int	adp,
		int	OPT,
		FILE *	sout
	)			`[private]`

double libMesh::VariationalMeshSmoother::localP	(	int	n,
		LPLPLPDOUBLE	W,
		LPLPDOUBLE	F,
		LPLPDOUBLE	R,
		LPINT	cell_in,
		LPINT	mask,
		double	epsilon,
		double	w,
		int	nvert,
		LPLPDOUBLE	H,
		int	me,
		double	vol,
		int	f,
		double *	Vmin,
		double *	qmin,
		int	adp,
		LPDOUBLE	afun,
		LPDOUBLE	Gloc,
		FILE *	sout
	)			`[private]`

composes local matrix W and right side F from all quadrature nodes of one cell

Definition at line 1952 of file mesh_smoother_vsmoother.C.

References avertex(), and vertex().

Referenced by minJ(), and minJ_BC().

01955 {
01956 
01957   int  ii, jj, kk, i, j, k, l, m, nn;
01958   double  sigma=0.0, fun, lqmin, gqmin, g;
01959   double K[9];
01960   /*f - flag, f=0 for determination of Hessian and gradient of the functional,
01961     f=1 for determination of the functional value only;
01962     K - determines approximation rule for local integral over the cell*/
01963 
01964 
01965   /*-------adaptivity, determined on the first step for adp>0 (nodal based)------*/
01966   if(f==0)
01967   {
01968     if(adp>0)
01969       avertex(n, afun, Gloc, R, cell_in, nvert, adp, sout);
01970     if(adp==0)
01971     {
01972       for(i=0;i<n;i++)
01973         Gloc[i]=1.0;
01974     }
01975   }
01976 
01977 fun=0; gqmin=1e32; g=0;//Vmin
01978 
01979 //cell integration depending on cell type
01980 if(n==2){//2D
01981         if(nvert==3){//tri
01982                 sigma=1.0;
01983         fun+=vertex(n, W, F, R, cell_in, epsilon, w, nvert, K, H, me, vol, f, &lqmin, adp, Gloc, sigma, sout);
01984                 g+=sigma*lqmin;
01985         if(gqmin>lqmin) gqmin=lqmin;
01986         }
01987     if(nvert==4){//quad
01988         for(i=0; i<2; i++){ K[0]=i;
01989             for(j=0; j<2; j++){ K[1]=j;
01990                             sigma=0.25;
01991                             fun+=vertex(n, W, F, R, cell_in, epsilon, w, nvert, K, H, me,
01992                                 vol, f, &lqmin, adp, Gloc, sigma, sout);
01993                 g+=sigma*lqmin;
01994                 if(gqmin>lqmin) gqmin=lqmin;
01995             }
01996           }
01997         } else {//quad tri
01998             for(i=0; i<3; i++){ K[0]=i*0.5; k=i/2; K[1]=(double)k;
01999                for(j=0; j<3; j++){  K[2]=j*0.5; k=j/2; K[3]=(double)k;
02000                                if(i==j) sigma=1.0/12; else sigma=1.0/24;
02001                                fun+=vertex(n, W, F, R, cell_in, epsilon, w, nvert, K, H, me,
02002                                    vol, f, &lqmin, adp, Gloc, sigma, sout);
02003                            g+=sigma*lqmin;
02004                    if(gqmin>lqmin) gqmin=lqmin;
02005                            }
02006                         }
02007         }
02008 }
02009 if(n==3){//3D
02010        if(nvert==4){//tetr
02011                    sigma=1.0;
02012                    fun+=vertex(n, W, F, R, cell_in, epsilon, w, nvert, K, H, me,
02013                            vol, f, &lqmin, adp, Gloc, sigma, sout);
02014                    g+=sigma*lqmin;
02015            if(gqmin>lqmin) gqmin=lqmin;
02016        }
02017        if(nvert==6){//prism
02018           for(i=0;i<2;i++){ K[0]=i;
02019               for(j=0;j<2;j++){ K[1]=j;
02020                   for(k=0;k<3;k++) {K[2]=(double)k/2.0; K[3]=(double)(k%2);
02021                                       sigma=1.0/12.0;
02022                                       fun+=vertex(n, W, F, R, cell_in, epsilon, w, nvert, K, H, me,
02023                                       vol, f, &lqmin, adp, Gloc, sigma, sout);
02024                               g+=sigma*lqmin;
02025                       if(gqmin>lqmin) gqmin=lqmin;
02026                                   }
02027                           }
02028                   }
02029            }
02030        if(nvert==8){//hex
02031          for(i=0; i<2; i++){ K[0]=i;
02032              for(j=0; j<2; j++){ K[1]=j;
02033                  for(k=0; k<2; k++){ K[2]=k;
02034                      for(l=0; l<2; l++){ K[3]=l;
02035                          for(m=0; m<2; m++){ K[4]=m;
02036                              for(nn=0; nn<2; nn++){ K[5]=nn;
02037                                                             if((i==nn)&&(j==l)&&(k==m)) sigma=(double)1/27;
02038                                 if(((i==nn)&&(j==l)&&(k!=m))||((i==nn)&&(j!=l)&&(k==m))||((i!=nn)&&(j==l)&&(k==m))) sigma=(double)1/(27*2);
02039                                 if(((i==nn)&&(j!=l)&&(k!=m))||((i!=nn)&&(j!=l)&&(k==m))||((i!=nn)&&(j==l)&&(k!=m))) sigma=(double)1/(27*4);
02040                                 if((i!=nn)&&(j!=l)&&(k!=m)) sigma=(double)1/(27*8);
02041                                                             fun+=vertex(n, W, F, R, cell_in, epsilon, w, nvert, K, H, me,
02042                                                 vol, f, &lqmin, adp, Gloc, sigma, sout);
02043                                         g+=sigma*lqmin;
02044                                 if(gqmin>lqmin) gqmin=lqmin;
02045                                                          }
02046                                                  }
02047                                          }
02048                                  }
02049                          }
02050                  }
02051            } else {//quad tetr
02052           for(i=0;i<4;i++){
02053             for(j=0;j<4;j++){
02054                   for(k=0;k<4;k++){
02055                     switch (i)
02056                     { case 0 : K[0]=0; K[1]=0; K[2]=0; break;
02057                       case 1 : K[0]=1; K[1]=0; K[2]=0; break;
02058                       case 2 : K[0]=1.0/2; K[1]=1; K[2]=0; break;
02059                       case 3 : K[0]=1.0/2; K[1]=1.0/3; K[2]=1; break; }
02060                     switch (j)
02061                     { case 0 : K[3]=0; K[4]=0; K[5]=0; break;
02062                       case 1 : K[3]=1; K[4]=0; K[5]=0; break;
02063                       case 2 : K[3]=1.0/2; K[4]=1; K[5]=0; break;
02064                       case 3 : K[3]=1.0/2; K[4]=1.0/3; K[5]=1; break; }
02065                     switch (k)
02066                     { case 0 : K[6]=0; K[7]=0; K[8]=0; break;
02067                       case 1 : K[6]=1; K[7]=0; K[8]=0; break;
02068                       case 2 : K[6]=1.0/2; K[7]=1; K[8]=0; break;
02069                       case 3 : K[6]=1.0/2; K[7]=1.0/3; K[8]=1; break; }
02070                            if((i==j)&&(j==k)) sigma=1.0/120; else
02071                        if((i==j)||(j==k)||(i==k)) sigma=1.0/360; else sigma=1.0/720;
02072                            fun+=vertex(n, W, F, R, cell_in, epsilon, w, nvert, K, H, me,
02073                                vol, f, &lqmin, adp, Gloc, sigma, sout);
02074                        g+=sigma*lqmin;
02075                if(gqmin>lqmin) gqmin=lqmin;
02076                           }
02077                         }
02078                   }
02079            }
02080 }
02081 
02082 /*---fixed nodes correction---*/
02083 for(ii=0;ii<nvert;ii++)
02084 {
02085     if(mask[cell_in[ii]]==1)
02086     {
02087       for(kk=0;kk<n;kk++)
02088       {
02089         for(jj=0;jj<nvert;jj++)
02090         {
02091           W[kk][ii][jj]=0;
02092           W[kk][jj][ii]=0;
02093         }
02094 
02095         W[kk][ii][ii]=1;
02096         F[kk][ii]=0;
02097       }
02098     }
02099 }
02100 /*---end of fixed nodes correction---*/
02101 
02102 (*Vmin)=g;
02103 (*qmin)=gqmin/vol;
02104 
02105 return fun;
02106 }

double libMesh::VariationalMeshSmoother::maxE	(	int	n,
		int	N,
		LPLPDOUBLE	R,
		int	ncells,
		LPLPINT	cells,
		LPINT	mcells,
		int	me,
		LPLPLPDOUBLE	H,
		double	v,
		double	epsilon,
		double	w,
		LPDOUBLE	Gamma,
		double *	qmin,
		FILE *	sout
	)			`[private]`

Determines the values of maxE_theta

Definition at line 959 of file mesh_smoother_vsmoother.C.

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

00962 {
00963   LPLPDOUBLE Q;
00964   double K[9], a1[3], a2[3], a3[3];
00965   double  gemax, det, tr, E=0.0, sigma=0.0, chi, vmin;
00966   int  ii, i, j, k, l, m, nn, kk, ll;
00967 
00968   Q = alloc_d_n1_n2(3, 10);
00969   for(i=0;i<n;i++) Q[i]=alloc_d_n1(3*n+n%2);
00970 
00971   gemax=-1e32; vmin=1e32;
00972 
00973 for(ii=0; ii<ncells; ii++) if(mcells[ii]>=0){
00974     if(n==2){
00975                 if(cells[ii][3]==-1){//tri
00976           basisA(2,Q,3,K,H[ii],me);
00977           for(k=0;k<2;k++){a1[k]=0; a2[k]=0;
00978              for(l=0;l<3;l++){
00979                  a1[k]=a1[k]+Q[k][l]*R[cells[ii][l]][0];
00980                  a2[k]=a2[k]+Q[k][l]*R[cells[ii][l]][1];
00981              }
00982           }
00983           det=jac2(a1[0],a1[1],a2[0],a2[1]);
00984           tr=0.5*(a1[0]*a1[0]+a2[0]*a2[0]+a1[1]*a1[1]+a2[1]*a2[1]);
00985           chi=0.5*(det+sqrt(det*det+epsilon*epsilon));
00986           E=(1-w)*tr/chi+0.5*w*(v+det*det/v)/chi;
00987           if(E>gemax) gemax=E;
00988                   if(vmin>det) vmin=det;
00989                 } else if(cells[ii][4]==-1){ E=0; //quad
00990           for(i=0; i<2; i++){ K[0]=i;
00991             for(j=0; j<2; j++){ K[1]=j;
00992                            basisA(2,Q,4,K,H[ii],me);
00993                for(k=0;k<2;k++){a1[k]=0; a2[k]=0;
00994                   for(l=0;l<4;l++){
00995                       a1[k]=a1[k]+Q[k][l]*R[cells[ii][l]][0];
00996                       a2[k]=a2[k]+Q[k][l]*R[cells[ii][l]][1];
00997                   }
00998                }
00999                det=jac2(a1[0],a1[1],a2[0],a2[1]);
01000                tr=0.5*(a1[0]*a1[0]+a2[0]*a2[0]+a1[1]*a1[1]+a2[1]*a2[1]);
01001                chi=0.5*(det+sqrt(det*det+epsilon*epsilon));
01002                E+=0.25*((1-w)*tr/chi+0.5*w*(v+det*det/v)/chi);
01003                            if(vmin>det) vmin=det;
01004                         }
01005           }
01006           if(E>gemax) gemax=E;
01007                 } else { E=0; //quad tri
01008             for(i=0; i<3; i++){ K[0]=i*0.5; k=i/2; K[1]=(double)k;
01009                for(j=0; j<3; j++){  K[2]=j*0.5; k=j/2; K[3]=(double)k;
01010                               basisA(2,Q,6,K,H[ii],me);
01011                   for(k=0;k<2;k++){a1[k]=0; a2[k]=0;
01012                      for(l=0;l<6;l++){
01013                         a1[k]=a1[k]+Q[k][l]*R[cells[ii][l]][0];
01014                         a2[k]=a2[k]+Q[k][l]*R[cells[ii][l]][1];
01015                      }
01016                   }
01017                   det=jac2(a1[0],a1[1],a2[0],a2[1]);
01018                   if(i==j) sigma=1.0/12; else sigma=1.0/24;
01019                   tr=0.5*(a1[0]*a1[0]+a2[0]*a2[0]+a1[1]*a1[1]+a2[1]*a2[1]);
01020                   chi=0.5*(det+sqrt(det*det+epsilon*epsilon));
01021                   E+=sigma*((1-w)*tr/chi+0.5*w*(v+det*det/v)/chi);
01022                                   if(vmin>det) vmin=det;
01023                            }
01024                         }
01025             if(E>gemax) gemax=E;
01026                 }
01027         }
01028     if(n==3){
01029        if(cells[ii][4]==-1){//tetr
01030            basisA(3,Q,4,K,H[ii],me);
01031            for(k=0;k<3;k++){a1[k]=0; a2[k]=0; a3[k]=0;
01032               for(l=0;l<4;l++){
01033                   a1[k]=a1[k]+Q[k][l]*R[cells[ii][l]][0];
01034                   a2[k]=a2[k]+Q[k][l]*R[cells[ii][l]][1];
01035                   a3[k]=a3[k]+Q[k][l]*R[cells[ii][l]][2];
01036               }
01037            }
01038            det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
01039            tr=0; for(k=0;k<3;k++) tr+=(a1[k]*a1[k]+a2[k]*a2[k]+a3[k]*a3[k])/3.0;
01040            chi=0.5*(det+sqrt(det*det+epsilon*epsilon));
01041            E=(1-w)*pow(tr,1.5)/chi+0.5*w*(v+det*det/v)/chi;
01042            if(E>gemax) gemax=E;
01043                    if(vmin>det) vmin=det;
01044        } else if(cells[ii][6]==-1){//prism
01045           E=0;
01046           for(i=0;i<2;i++){ K[0]=i;
01047               for(j=0;j<2;j++){ K[1]=j;
01048                   for(k=0;k<3;k++) {K[2]=(double)k/2.0; K[3]=(double)(k%2);
01049                       basisA(3,Q,6,K,H[ii],me);
01050                       for(kk=0;kk<3;kk++){a1[kk]=0; a2[kk]=0; a3[kk]=0;
01051                           for(ll=0;ll<6;ll++){
01052                               a1[kk]=a1[kk]+Q[kk][ll]*R[cells[ii][ll]][0];
01053                               a2[kk]=a2[kk]+Q[kk][ll]*R[cells[ii][ll]][1];
01054                               a3[kk]=a3[kk]+Q[kk][ll]*R[cells[ii][ll]][2];
01055                           }
01056                       }
01057                       det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
01058                       tr=0; for(kk=0;kk<3;kk++) tr+=(a1[kk]*a1[kk]+a2[kk]*a2[kk]+a3[kk]*a3[kk])/3.0;
01059                       chi=0.5*(det+sqrt(det*det+epsilon*epsilon));
01060                       E+=((1-w)*pow(tr,1.5)/chi+0.5*w*(v+det*det/v)/chi)/12.0;
01061                                           if(vmin>det) vmin=det;
01062                                   }
01063                           }
01064                   }
01065           if(E>gemax) gemax=E;
01066        } else if(cells[ii][8]==-1){ E=0; //hex
01067          for(i=0; i<2; i++){ K[0]=i;
01068              for(j=0; j<2; j++){ K[1]=j;
01069                  for(k=0; k<2; k++){ K[2]=k;
01070                      for(l=0; l<2; l++){ K[3]=l;
01071                          for(m=0; m<2; m++){ K[4]=m;
01072                              for(nn=0; nn<2; nn++){ K[5]=nn;
01073                                 basisA(3,Q,8,K,H[ii],me);
01074                                 for(kk=0;kk<3;kk++){a1[kk]=0; a2[kk]=0; a3[kk]=0;
01075                                    for(ll=0;ll<8;ll++){
01076                                       a1[kk]=a1[kk]+Q[kk][ll]*R[cells[ii][ll]][0];
01077                                       a2[kk]=a2[kk]+Q[kk][ll]*R[cells[ii][ll]][1];
01078                                       a3[kk]=a3[kk]+Q[kk][ll]*R[cells[ii][ll]][2];
01079                                    }
01080                                 }
01081                                 det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
01082                                                                 if((i==nn)&&(j==l)&&(k==m)) sigma=(double)1/27;
01083                                 if(((i==nn)&&(j==l)&&(k!=m))||((i==nn)&&(j!=l)&&(k==m))||((i!=nn)&&(j==l)&&(k==m))) sigma=(double)1/(27*2);
01084                                 if(((i==nn)&&(j!=l)&&(k!=m))||((i!=nn)&&(j!=l)&&(k==m))||((i!=nn)&&(j==l)&&(k!=m))) sigma=(double)1/(27*4);
01085                                 if((i!=nn)&&(j!=l)&&(k!=m)) sigma=(double)1/(27*8);
01086                                 tr=0; for(kk=0;kk<3;kk++) tr+=(a1[kk]*a1[kk]+a2[kk]*a2[kk]+a3[kk]*a3[kk])/3.0;
01087                                 chi=0.5*(det+sqrt(det*det+epsilon*epsilon));
01088                                 E+=((1-w)*pow(tr,1.5)/chi+0.5*w*(v+det*det/v)/chi)*sigma;
01089                                                                 if(vmin>det) vmin=det;
01090                                                          }
01091                          }
01092                       }
01093                   }
01094               }
01095           }
01096           if(E>gemax) gemax=E;
01097        } else {//quad tetr
01098                    E=0;
01099           for(i=0;i<4;i++){
01100             for(j=0;j<4;j++){
01101                   for(k=0;k<4;k++){
01102                     switch (i)
01103                     { case 0 : K[0]=0; K[1]=0; K[2]=0; break;
01104                       case 1 : K[0]=1; K[1]=0; K[2]=0; break;
01105                       case 2 : K[0]=1.0/2; K[1]=1; K[2]=0; break;
01106                       case 3 : K[0]=1.0/2; K[1]=1.0/3; K[2]=1; break; }
01107                     switch (j)
01108                     { case 0 : K[3]=0; K[4]=0; K[5]=0; break;
01109                       case 1 : K[3]=1; K[4]=0; K[5]=0; break;
01110                       case 2 : K[3]=1.0/2; K[4]=1; K[5]=0; break;
01111                       case 3 : K[3]=1.0/2; K[4]=1.0/3; K[5]=1; break; }
01112                     switch (k)
01113                     { case 0 : K[6]=0; K[7]=0; K[8]=0; break;
01114                       case 1 : K[6]=1; K[7]=0; K[8]=0; break;
01115                       case 2 : K[6]=1.0/2; K[7]=1; K[8]=0; break;
01116                       case 3 : K[6]=1.0/2; K[7]=1.0/3; K[8]=1; break; }
01117                     basisA(3,Q,10,K,H[ii],me);
01118                     for(kk=0;kk<3;kk++){a1[kk]=0; a2[kk]=0; a3[kk]=0;
01119                        for(ll=0;ll<10;ll++){
01120                           a1[kk]=a1[kk]+Q[kk][ll]*R[cells[ii][ll]][0];
01121                           a2[kk]=a2[kk]+Q[kk][ll]*R[cells[ii][ll]][1];
01122                           a3[kk]=a3[kk]+Q[kk][ll]*R[cells[ii][ll]][2];
01123                        }
01124                     }
01125                     det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
01126                     if((i==j)&&(j==k)) sigma=1.0/120; else
01127                             if((i==j)||(j==k)||(i==k)) sigma=1.0/360; else sigma=1.0/720;
01128                     tr=0; for(kk=0;kk<3;kk++) tr+=(a1[kk]*a1[kk]+a2[kk]*a2[kk]+a3[kk]*a3[kk])/3.0;
01129                     chi=0.5*(det+sqrt(det*det+epsilon*epsilon));
01130                     E+=((1-w)*pow(tr,1.5)/chi+0.5*w*(v+det*det/v)/chi)*sigma;
01131                                         if(vmin>det) vmin=det;
01132                           }
01133                         }
01134                   }
01135           if(E>gemax) gemax=E;
01136            }
01137         }
01138         Gamma[ii]=E;
01139 }
01140 
01141 (*qmin)=vmin;
01142 
01143 for(i=0;i<n;i++) free(Q[i]);
01144 free(Q);
01145 
01146 return gemax;
01147 }

void libMesh::VariationalMeshSmoother::metr_data_gen	(	char	grid[],
		char	metr[],
		int	n,
		int	me,
		FILE *	sout
	)			`[private]`

Metric Generration

Definition at line 2627 of file mesh_smoother_vsmoother.C.

References alloc_d_n1(), alloc_d_n1_n2(), alloc_i_n1(), alloc_i_n1_n2(), basisA(), jac2(), jac3(), std::pow(), and readgr().

Referenced by smooth().

02628 {
02629   LPLPDOUBLE R;
02630   LPLPINT cells;
02631   LPINT mask, mcells;
02632   LPLPDOUBLE Q;
02633   double K[9], a1[3], a2[3], a3[3];
02634   double det, g1, g2, g3, det_o, g1_o, g2_o, g3_o, eps=1e-3;
02635   int  i, j, k, l, N, ncells, Ncells, nvert, scanned;
02636   FILE *stream;
02637 
02638   Q = alloc_d_n1_n2(3, 10);
02639   for(i=0;i<n;i++) Q[i]=alloc_d_n1(3*n+n%2);
02640 
02641   //read the initial mesh
02642   stream=fopen(grid,"r");
02643   scanned = fscanf(stream, "%d \n%d \n%d \n%d \n",&i,&N,&ncells,&j);
02644   libmesh_assert_not_equal_to (scanned, EOF);
02645   fclose(stream);
02646 
02647   mask=alloc_i_n1(N);
02648   mcells=alloc_i_n1(ncells);
02649   R=alloc_d_n1_n2(N,n);
02650   cells=alloc_i_n1_n2(ncells,3*n+n%2);
02651   for(i=0;i<ncells;i++) cells[i]=alloc_i_n1(3*n+n%2);
02652   for(i=0;i<N;i++) R[i]=alloc_d_n1(n);
02653 
02654   readgr(n,N,R,mask,ncells,cells,mcells,0,mcells,mcells, sout);
02655 
02656   //genetrate metric file
02657   stream=fopen(metr,"w+");
02658   Ncells=0;
02659   det_o=1.0; g1_o=1.0; g2_o=1.0; g3_o=1.0;
02660   for(i=0; i<ncells; i++) if(mcells[i]>=0){
02661           nvert=0; while(cells[i][nvert]>=0) nvert++;
02662           if(n==2){//2D - tri and quad
02663                   if(nvert==3){//tri
02664                           basisA(2,Q,3,K,Q,1);
02665                           for(k=0;k<2;k++){a1[k]=0; a2[k]=0;
02666                               for(l=0;l<3;l++){
02667                                       a1[k]=a1[k]+Q[k][l]*R[cells[i][l]][0];
02668                                       a2[k]=a2[k]+Q[k][l]*R[cells[i][l]][1];
02669                                   }
02670                           }
02671                           det=jac2(a1[0],a1[1],a2[0],a2[1]);
02672                           g1=sqrt(a1[0]*a1[0]+a2[0]*a2[0]);
02673                           g2=sqrt(a1[1]*a1[1]+a2[1]*a2[1]);
02674                           //need to keep data from previous cell
02675                           if((fabs(det)<eps*eps*det_o)||(det<0)) det=det_o;
02676                           if((fabs(g1)<eps*g1_o)||(g1<0)) g1=g1_o;
02677                           if((fabs(g2)<eps*g2_o)||(g2<0)) g2=g2_o;
02678                           //write to file
02679                           if(me==2) fprintf(stream,"%e 0.000000e+00 \n0.000000e+00 %e \n",1.0/sqrt(det),1.0/sqrt(det));
02680               if(me==3) fprintf(stream,"%e 0.000000e+00 \n0.000000e+00 %e \n",1.0/g1,1.0/g2);
02681                           det_o=det; g1_o=g1; g2_o=g2;
02682                           Ncells++;
02683                   }
02684                   if(nvert==4){//quad
02685               a1[0]=R[cells[i][1]][0]-R[cells[i][0]][0];
02686               a1[1]=R[cells[i][2]][0]-R[cells[i][0]][0];
02687               a2[0]=R[cells[i][1]][1]-R[cells[i][0]][1];
02688               a2[1]=R[cells[i][2]][1]-R[cells[i][0]][1];
02689                           det=jac2(a1[0],a1[1],a2[0],a2[1]);
02690                           g1=sqrt(a1[0]*a1[0]+a2[0]*a2[0]);
02691                           g2=sqrt(a1[1]*a1[1]+a2[1]*a2[1]);
02692                           //need to keep data from previous cell
02693                           if((fabs(det)<eps*eps*det_o)||(det<0)) det=det_o;
02694                           if((fabs(g1)<eps*g1_o)||(g1<0)) g1=g1_o;
02695                           if((fabs(g2)<eps*g2_o)||(g2<0)) g2=g2_o;
02696                           //write to file
02697                           if(me==2) fprintf(stream,"%e %e \n0.000000e+00 %e \n",1.0/sqrt(det),0.5/sqrt(det),0.5*sqrt(3.0/det));
02698               if(me==3) fprintf(stream,"%e %e \n0.000000e+00 %e \n",1.0/g1,0.5/g2,0.5*sqrt(3.0)/g2);
02699                           det_o=det; g1_o=g1; g2_o=g2;
02700                           Ncells++;
02701                       a1[0]=R[cells[i][3]][0]-R[cells[i][1]][0];
02702                       a1[1]=R[cells[i][0]][0]-R[cells[i][1]][0];
02703                       a2[0]=R[cells[i][3]][1]-R[cells[i][1]][1];
02704                       a2[1]=R[cells[i][0]][1]-R[cells[i][1]][1];
02705                           det=jac2(a1[0],a1[1],a2[0],a2[1]);
02706                           g1=sqrt(a1[0]*a1[0]+a2[0]*a2[0]);
02707                           g2=sqrt(a1[1]*a1[1]+a2[1]*a2[1]);
02708                           //need to keep data from previous cell
02709                           if((fabs(det)<eps*eps*det_o)||(det<0)) det=det_o;
02710                           if((fabs(g1)<eps*g1_o)||(g1<0)) g1=g1_o;
02711                           if((fabs(g2)<eps*g2_o)||(g2<0)) g2=g2_o;
02712                           //write to file
02713                           if(me==2) fprintf(stream,"%e %e \n0.000000e+00 %e \n",1.0/sqrt(det),0.5/sqrt(det),0.5*sqrt(3.0/det));
02714               if(me==3) fprintf(stream,"%e %e \n0.000000e+00 %e \n",1.0/g1,0.5/g2,0.5*sqrt(3.0)/g2);
02715                           det_o=det; g1_o=g1; g2_o=g2;
02716                           Ncells++;
02717                       a1[0]=R[cells[i][2]][0]-R[cells[i][3]][0];
02718                       a1[1]=R[cells[i][1]][0]-R[cells[i][3]][0];
02719                       a2[0]=R[cells[i][2]][1]-R[cells[i][3]][1];
02720                       a2[1]=R[cells[i][1]][1]-R[cells[i][3]][1];
02721                           det=jac2(a1[0],a1[1],a2[0],a2[1]);
02722                           g1=sqrt(a1[0]*a1[0]+a2[0]*a2[0]);
02723                           g2=sqrt(a1[1]*a1[1]+a2[1]*a2[1]);
02724                           //need to keep data from previous cell
02725                           if((fabs(det)<eps*eps*det_o)||(det<0)) det=det_o;
02726                           if((fabs(g1)<eps*g1_o)||(g1<0)) g1=g1_o;
02727                           if((fabs(g2)<eps*g2_o)||(g2<0)) g2=g2_o;
02728                           //write to file
02729                           if(me==2) fprintf(stream,"%e %e \n0.000000e+00 %e \n",1.0/sqrt(det),0.5/sqrt(det),0.5*sqrt(3.0/det));
02730               if(me==3) fprintf(stream,"%e %e \n0.000000e+00 %e \n",1.0/g1,0.5/g2,0.5*sqrt(3.0)/g2);
02731                           det_o=det; g1_o=g1; g2_o=g2;
02732                           Ncells++;
02733                       a1[0]=R[cells[i][0]][0]-R[cells[i][2]][0];
02734                       a1[1]=R[cells[i][3]][0]-R[cells[i][2]][0];
02735                       a2[0]=R[cells[i][0]][1]-R[cells[i][2]][1];
02736                       a2[1]=R[cells[i][3]][1]-R[cells[i][2]][1];
02737                           det=jac2(a1[0],a1[1],a2[0],a2[1]);
02738                           g1=sqrt(a1[0]*a1[0]+a2[0]*a2[0]);
02739                           g2=sqrt(a1[1]*a1[1]+a2[1]*a2[1]);
02740                           //need to keep data from previous cell
02741                           if((fabs(det)<eps*eps*det_o)||(det<0)) det=det_o;
02742                           if((fabs(g1)<eps*g1_o)||(g1<0)) g1=g1_o;
02743                           if((fabs(g2)<eps*g2_o)||(g2<0)) g2=g2_o;
02744                           //write to file
02745                           if(me==2) fprintf(stream,"%e %e \n0.000000e+00 %e \n",1.0/sqrt(det),0.5/sqrt(det),0.5*sqrt(3.0/det));
02746               if(me==3) fprintf(stream,"%e %e \n0.000000e+00 %e \n",1.0/g1,0.5/g2,0.5*sqrt(3.0)/g2);
02747                           det_o=det; g1_o=g1; g2_o=g2;
02748                           Ncells++;
02749 
02750                   }
02751           }
02752           if(n==3){//3D - tetr and hex
02753                   if(nvert==4){//tetr
02754                           basisA(3,Q,4,K,Q,1);
02755                           for(k=0;k<3;k++){a1[k]=0; a2[k]=0; a3[k]=0;
02756                               for(l=0;l<4;l++){
02757                                           a1[k]=a1[k]+Q[k][l]*R[cells[i][l]][0];
02758                                           a2[k]=a2[k]+Q[k][l]*R[cells[i][l]][1];
02759                                           a3[k]=a3[k]+Q[k][l]*R[cells[i][l]][2];
02760                                   }
02761                           }
02762                           det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
02763               g1=sqrt(a1[0]*a1[0]+a2[0]*a2[0]+a3[0]*a3[0]);
02764                           g2=sqrt(a1[1]*a1[1]+a2[1]*a2[1]+a3[1]*a3[1]);
02765                           g3=sqrt(a1[2]*a1[2]+a2[2]*a2[2]+a3[2]*a3[2]);
02766                           //need to keep data from previous cell
02767                           if((fabs(det)<eps*eps*eps*det_o)||(det<0)) det=det_o;
02768                           if((fabs(g1)<eps*g1_o)||(g1<0)) g1=g1_o;
02769                           if((fabs(g2)<eps*g2_o)||(g2<0)) g2=g2_o;
02770                           if((fabs(g3)<eps*g3_o)||(g3<0)) g3=g3_o;
02771                           //write to file
02772                           if(me==2) fprintf(stream,"%e 0.000000e+00  0.000000e+00\n0.000000e+00 %e 0.000000e+00\n  0.000000e+00 0.000000e+00 %e\n",1.0/pow(det,1.0/3.0),1.0/pow(det,1.0/3.0),1.0/pow(det,1.0/3.0));
02773               if(me==3) fprintf(stream,"%e 0.000000e+00  0.000000e+00\n0.000000e+00 %e 0.000000e+00\n  0.000000e+00 0.000000e+00 %e\n",1.0/g1,1.0/g2,1.0/g3);
02774                           det_o=det; g1_o=g1; g2_o=g2; g3_o=g3;
02775                           Ncells++;
02776                   }
02777                   if(nvert==8){//hex
02778                      a1[0]=R[cells[i][1]][0]-R[cells[i][0]][0];
02779                      a1[1]=R[cells[i][2]][0]-R[cells[i][0]][0];
02780                      a1[2]=R[cells[i][4]][0]-R[cells[i][0]][0];
02781                      a2[0]=R[cells[i][1]][1]-R[cells[i][0]][1];
02782                      a2[1]=R[cells[i][2]][1]-R[cells[i][0]][1];
02783                      a2[2]=R[cells[i][4]][1]-R[cells[i][0]][1];
02784                      a3[0]=R[cells[i][1]][2]-R[cells[i][0]][2];
02785                      a3[1]=R[cells[i][2]][2]-R[cells[i][0]][2];
02786                      a3[2]=R[cells[i][4]][2]-R[cells[i][0]][2];
02787                      det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
02788                      g1=sqrt(a1[0]*a1[0]+a2[0]*a2[0]+a3[0]*a3[0]);
02789                           g2=sqrt(a1[1]*a1[1]+a2[1]*a2[1]+a3[1]*a3[1]);
02790                           g3=sqrt(a1[2]*a1[2]+a2[2]*a2[2]+a3[2]*a3[2]);
02791                           //need to keep data from previous cell
02792                           if((fabs(det)<eps*eps*eps*det_o)||(det<0)) det=det_o;
02793                           if((fabs(g1)<eps*g1_o)||(g1<0)) g1=g1_o;
02794                           if((fabs(g2)<eps*g2_o)||(g2<0)) g2=g2_o;
02795                           if((fabs(g3)<eps*g3_o)||(g3<0)) g3=g3_o;
02796                           //write to file
02797                           if(me==2) fprintf(stream,"%e %e %e\n0.000000e+00 %e %e\n  0.000000e+00 0.000000e+00 %e\n",1.0/pow(det,1.0/3.0),0.5/pow(det,1.0/3.0),0.5/pow(det,1.0/3.0),0.5*sqrt(3.0)/pow(det,1.0/3.0),0.5/(sqrt(3.0)*pow(det,1.0/3.0)),sqrt(2/3.0)/pow(det,1.0/3.0));
02798                           if(me==3) fprintf(stream,"%e %e %e\n0.000000e+00 %e %e\n  0.000000e+00 0.000000e+00 %e\n",1.0/g1,0.5/g2,0.5/g3,0.5*sqrt(3.0)/g2,0.5/(sqrt(3.0)*g3),sqrt(2/3.0)/g3);
02799                           det_o=det; g1_o=g1; g2_o=g2; g3_o=g3;
02800                           Ncells++;
02801                      a1[0]=R[cells[i][3]][0]-R[cells[i][1]][0];
02802                      a1[1]=R[cells[i][0]][0]-R[cells[i][1]][0];
02803                      a1[2]=R[cells[i][5]][0]-R[cells[i][1]][0];
02804                      a2[0]=R[cells[i][3]][1]-R[cells[i][1]][1];
02805                      a2[1]=R[cells[i][0]][1]-R[cells[i][1]][1];
02806                      a2[2]=R[cells[i][5]][1]-R[cells[i][1]][1];
02807                      a3[0]=R[cells[i][3]][2]-R[cells[i][1]][2];
02808                      a3[1]=R[cells[i][0]][2]-R[cells[i][1]][2];
02809                      a3[2]=R[cells[i][5]][2]-R[cells[i][1]][2];
02810                      det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
02811                      g1=sqrt(a1[0]*a1[0]+a2[0]*a2[0]+a3[0]*a3[0]);
02812                           g2=sqrt(a1[1]*a1[1]+a2[1]*a2[1]+a3[1]*a3[1]);
02813                           g3=sqrt(a1[2]*a1[2]+a2[2]*a2[2]+a3[2]*a3[2]);
02814                           //need to keep data from previous cell
02815                           if((fabs(det)<eps*eps*eps*det_o)||(det<0)) det=det_o;
02816                           if((fabs(g1)<eps*g1_o)||(g1<0)) g1=g1_o;
02817                           if((fabs(g2)<eps*g2_o)||(g2<0)) g2=g2_o;
02818                           if((fabs(g3)<eps*g3_o)||(g3<0)) g3=g3_o;
02819                           //write to file
02820                           if(me==2) fprintf(stream,"%e %e %e\n0.000000e+00 %e %e\n  0.000000e+00 0.000000e+00 %e\n",1.0/pow(det,1.0/3.0),0.5/pow(det,1.0/3.0),0.5/pow(det,1.0/3.0),0.5*sqrt(3.0)/pow(det,1.0/3.0),0.5/(sqrt(3.0)*pow(det,1.0/3.0)),sqrt(2/3.0)/pow(det,1.0/3.0));
02821                           if(me==3) fprintf(stream,"%e %e %e\n0.000000e+00 %e %e\n  0.000000e+00 0.000000e+00 %e\n",1.0/g1,0.5/g2,0.5/g3,0.5*sqrt(3.0)/g2,0.5/(sqrt(3.0)*g3),sqrt(2/3.0)/g3);
02822                           det_o=det; g1_o=g1; g2_o=g2; g3_o=g3;
02823                           Ncells++;
02824                      a1[0]=R[cells[i][2]][0]-R[cells[i][3]][0];
02825                      a1[1]=R[cells[i][1]][0]-R[cells[i][3]][0];
02826                      a1[2]=R[cells[i][7]][0]-R[cells[i][3]][0];
02827                      a2[0]=R[cells[i][2]][1]-R[cells[i][3]][1];
02828                      a2[1]=R[cells[i][1]][1]-R[cells[i][3]][1];
02829                      a2[2]=R[cells[i][7]][1]-R[cells[i][3]][1];
02830                      a3[0]=R[cells[i][2]][2]-R[cells[i][3]][2];
02831                      a3[1]=R[cells[i][1]][2]-R[cells[i][3]][2];
02832                      a3[2]=R[cells[i][7]][2]-R[cells[i][3]][2];
02833                      det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
02834                      g1=sqrt(a1[0]*a1[0]+a2[0]*a2[0]+a3[0]*a3[0]);
02835                           g2=sqrt(a1[1]*a1[1]+a2[1]*a2[1]+a3[1]*a3[1]);
02836                           g3=sqrt(a1[2]*a1[2]+a2[2]*a2[2]+a3[2]*a3[2]);
02837                           //need to keep data from previous cell
02838                           if((fabs(det)<eps*eps*eps*det_o)||(det<0)) det=det_o;
02839                           if((fabs(g1)<eps*g1_o)||(g1<0)) g1=g1_o;
02840                           if((fabs(g2)<eps*g2_o)||(g2<0)) g2=g2_o;
02841                           if((fabs(g3)<eps*g3_o)||(g3<0)) g3=g3_o;
02842                           //write to file
02843                           if(me==2) fprintf(stream,"%e %e %e\n0.000000e+00 %e %e\n  0.000000e+00 0.000000e+00 %e\n",1.0/pow(det,1.0/3.0),0.5/pow(det,1.0/3.0),0.5/pow(det,1.0/3.0),0.5*sqrt(3.0)/pow(det,1.0/3.0),0.5/(sqrt(3.0)*pow(det,1.0/3.0)),sqrt(2/3.0)/pow(det,1.0/3.0));
02844                           if(me==3) fprintf(stream,"%e %e %e\n0.000000e+00 %e %e\n  0.000000e+00 0.000000e+00 %e\n",1.0/g1,0.5/g2,0.5/g3,0.5*sqrt(3.0)/g2,0.5/(sqrt(3.0)*g3),sqrt(2/3.0)/g3);
02845                           det_o=det; g1_o=g1; g2_o=g2; g3_o=g3;
02846                           Ncells++;
02847                      a1[0]=R[cells[i][0]][0]-R[cells[i][2]][0];
02848                      a1[1]=R[cells[i][3]][0]-R[cells[i][2]][0];
02849                      a1[2]=R[cells[i][6]][0]-R[cells[i][2]][0];
02850                      a2[0]=R[cells[i][0]][1]-R[cells[i][2]][1];
02851                      a2[1]=R[cells[i][3]][1]-R[cells[i][2]][1];
02852                      a2[2]=R[cells[i][6]][1]-R[cells[i][2]][1];
02853                      a3[0]=R[cells[i][0]][2]-R[cells[i][2]][2];
02854                      a3[1]=R[cells[i][3]][2]-R[cells[i][2]][2];
02855                      a3[2]=R[cells[i][6]][2]-R[cells[i][2]][2];
02856                      det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
02857                      g1=sqrt(a1[0]*a1[0]+a2[0]*a2[0]+a3[0]*a3[0]);
02858                           g2=sqrt(a1[1]*a1[1]+a2[1]*a2[1]+a3[1]*a3[1]);
02859                           g3=sqrt(a1[2]*a1[2]+a2[2]*a2[2]+a3[2]*a3[2]);
02860                           //need to keep data from previous cell
02861                           if((fabs(det)<eps*eps*eps*det_o)||(det<0)) det=det_o;
02862                           if((fabs(g1)<eps*g1_o)||(g1<0)) g1=g1_o;
02863                           if((fabs(g2)<eps*g2_o)||(g2<0)) g2=g2_o;
02864                           if((fabs(g3)<eps*g3_o)||(g3<0)) g3=g3_o;
02865                           //write to file
02866                           if(me==2) fprintf(stream,"%e %e %e\n0.000000e+00 %e %e\n  0.000000e+00 0.000000e+00 %e\n",1.0/pow(det,1.0/3.0),0.5/pow(det,1.0/3.0),0.5/pow(det,1.0/3.0),0.5*sqrt(3.0)/pow(det,1.0/3.0),0.5/(sqrt(3.0)*pow(det,1.0/3.0)),sqrt(2/3.0)/pow(det,1.0/3.0));
02867                           if(me==3) fprintf(stream,"%e %e %e\n0.000000e+00 %e %e\n  0.000000e+00 0.000000e+00 %e\n",1.0/g1,0.5/g2,0.5/g3,0.5*sqrt(3.0)/g2,0.5/(sqrt(3.0)*g3),sqrt(2/3.0)/g3);
02868                           det_o=det; g1_o=g1; g2_o=g2; g3_o=g3;
02869                           Ncells++;
02870                      a1[0]=R[cells[i][6]][0]-R[cells[i][4]][0];
02871                      a1[1]=R[cells[i][5]][0]-R[cells[i][4]][0];
02872                      a1[2]=R[cells[i][0]][0]-R[cells[i][4]][0];
02873                      a2[0]=R[cells[i][6]][1]-R[cells[i][4]][1];
02874                      a2[1]=R[cells[i][5]][1]-R[cells[i][4]][1];
02875                      a2[2]=R[cells[i][0]][1]-R[cells[i][4]][1];
02876                      a3[0]=R[cells[i][6]][2]-R[cells[i][4]][2];
02877                      a3[1]=R[cells[i][5]][2]-R[cells[i][4]][2];
02878                      a3[2]=R[cells[i][0]][2]-R[cells[i][4]][2];
02879                      det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
02880                      g1=sqrt(a1[0]*a1[0]+a2[0]*a2[0]+a3[0]*a3[0]);
02881                           g2=sqrt(a1[1]*a1[1]+a2[1]*a2[1]+a3[1]*a3[1]);
02882                           g3=sqrt(a1[2]*a1[2]+a2[2]*a2[2]+a3[2]*a3[2]);
02883                           //need to keep data from previous cell
02884                           if((fabs(det)<eps*eps*eps*det_o)||(det<0)) det=det_o;
02885                           if((fabs(g1)<eps*g1_o)||(g1<0)) g1=g1_o;
02886                           if((fabs(g2)<eps*g2_o)||(g2<0)) g2=g2_o;
02887                           if((fabs(g3)<eps*g3_o)||(g3<0)) g3=g3_o;
02888                           //write to file
02889                           if(me==2) fprintf(stream,"%e %e %e\n0.000000e+00 %e %e\n  0.000000e+00 0.000000e+00 %e\n",1.0/pow(det,1.0/3.0),0.5/pow(det,1.0/3.0),0.5/pow(det,1.0/3.0),0.5*sqrt(3.0)/pow(det,1.0/3.0),0.5/(sqrt(3.0)*pow(det,1.0/3.0)),sqrt(2/3.0)/pow(det,1.0/3.0));
02890                           if(me==3) fprintf(stream,"%e %e %e\n0.000000e+00 %e %e\n  0.000000e+00 0.000000e+00 %e\n",1.0/g1,0.5/g2,0.5/g3,0.5*sqrt(3.0)/g2,0.5/(sqrt(3.0)*g3),sqrt(2/3.0)/g3);
02891                           det_o=det; g1_o=g1; g2_o=g2; g3_o=g3;
02892                           Ncells++;
02893                      a1[0]=R[cells[i][4]][0]-R[cells[i][5]][0];
02894                      a1[1]=R[cells[i][7]][0]-R[cells[i][5]][0];
02895                      a1[2]=R[cells[i][1]][0]-R[cells[i][5]][0];
02896                      a2[0]=R[cells[i][4]][1]-R[cells[i][5]][1];
02897                      a2[1]=R[cells[i][7]][1]-R[cells[i][5]][1];
02898                      a2[2]=R[cells[i][1]][1]-R[cells[i][5]][1];
02899                      a3[0]=R[cells[i][4]][2]-R[cells[i][5]][2];
02900                      a3[1]=R[cells[i][7]][2]-R[cells[i][5]][2];
02901                      a3[2]=R[cells[i][1]][2]-R[cells[i][5]][2];
02902                      det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
02903                      g1=sqrt(a1[0]*a1[0]+a2[0]*a2[0]+a3[0]*a3[0]);
02904                           g2=sqrt(a1[1]*a1[1]+a2[1]*a2[1]+a3[1]*a3[1]);
02905                           g3=sqrt(a1[2]*a1[2]+a2[2]*a2[2]+a3[2]*a3[2]);
02906                           //need to keep data from previous cell
02907                           if((fabs(det)<eps*eps*eps*det_o)||(det<0)) det=det_o;
02908                           if((fabs(g1)<eps*g1_o)||(g1<0)) g1=g1_o;
02909                           if((fabs(g2)<eps*g2_o)||(g2<0)) g2=g2_o;
02910                           if((fabs(g3)<eps*g3_o)||(g3<0)) g3=g3_o;
02911                           //write to file
02912                           if(me==2) fprintf(stream,"%e %e %e\n0.000000e+00 %e %e\n  0.000000e+00 0.000000e+00 %e\n",1.0/pow(det,1.0/3.0),0.5/pow(det,1.0/3.0),0.5/pow(det,1.0/3.0),0.5*sqrt(3.0)/pow(det,1.0/3.0),0.5/(sqrt(3.0)*pow(det,1.0/3.0)),sqrt(2/3.0)/pow(det,1.0/3.0));
02913                           if(me==3) fprintf(stream,"%e %e %e\n0.000000e+00 %e %e\n  0.000000e+00 0.000000e+00 %e\n",1.0/g1,0.5/g2,0.5/g3,0.5*sqrt(3.0)/g2,0.5/(sqrt(3.0)*g3),sqrt(2/3.0)/g3);
02914                           det_o=det; g1_o=g1; g2_o=g2; g3_o=g3;
02915                           Ncells++;
02916                      a1[0]=R[cells[i][5]][0]-R[cells[i][7]][0];
02917                      a1[1]=R[cells[i][6]][0]-R[cells[i][7]][0];
02918                      a1[2]=R[cells[i][3]][0]-R[cells[i][7]][0];
02919                      a2[0]=R[cells[i][5]][1]-R[cells[i][7]][1];
02920                      a2[1]=R[cells[i][6]][1]-R[cells[i][7]][1];
02921                      a2[2]=R[cells[i][3]][1]-R[cells[i][7]][1];
02922                      a3[0]=R[cells[i][5]][2]-R[cells[i][7]][2];
02923                      a3[1]=R[cells[i][6]][2]-R[cells[i][7]][2];
02924                      a3[2]=R[cells[i][3]][2]-R[cells[i][7]][2];
02925                      det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
02926                      g1=sqrt(a1[0]*a1[0]+a2[0]*a2[0]+a3[0]*a3[0]);
02927                           g2=sqrt(a1[1]*a1[1]+a2[1]*a2[1]+a3[1]*a3[1]);
02928                           g3=sqrt(a1[2]*a1[2]+a2[2]*a2[2]+a3[2]*a3[2]);
02929                           //need to keep data from previous cell
02930                           if((fabs(det)<eps*eps*eps*det_o)||(det<0)) det=det_o;
02931                           if((fabs(g1)<eps*g1_o)||(g1<0)) g1=g1_o;
02932                           if((fabs(g2)<eps*g2_o)||(g2<0)) g2=g2_o;
02933                           if((fabs(g3)<eps*g3_o)||(g3<0)) g3=g3_o;
02934                           //write to file
02935                           if(me==2) fprintf(stream,"%e %e %e\n0.000000e+00 %e %e\n  0.000000e+00 0.000000e+00 %e\n",1.0/pow(det,1.0/3.0),0.5/pow(det,1.0/3.0),0.5/pow(det,1.0/3.0),0.5*sqrt(3.0)/pow(det,1.0/3.0),0.5/(sqrt(3.0)*pow(det,1.0/3.0)),sqrt(2/3.0)/pow(det,1.0/3.0));
02936               if(me==3) fprintf(stream,"%e %e %e\n0.000000e+00 %e %e\n  0.000000e+00 0.000000e+00 %e\n",1.0/g1,0.5/g2,0.5/g3,0.5*sqrt(3.0)/g2,0.5/(sqrt(3.0)*g3),sqrt(2/3.0)/g3);
02937                           det_o=det; g1_o=g1; g2_o=g2; g3_o=g3;
02938                           Ncells++;
02939                      a1[0]=R[cells[i][6]][0]-R[cells[i][6]][0];
02940                      a1[1]=R[cells[i][4]][0]-R[cells[i][6]][0];
02941                      a1[2]=R[cells[i][2]][0]-R[cells[i][6]][0];
02942                      a2[0]=R[cells[i][6]][1]-R[cells[i][6]][1];
02943                      a2[1]=R[cells[i][4]][1]-R[cells[i][6]][1];
02944                      a2[2]=R[cells[i][2]][1]-R[cells[i][6]][1];
02945                      a3[0]=R[cells[i][6]][2]-R[cells[i][6]][2];
02946                      a3[1]=R[cells[i][4]][2]-R[cells[i][6]][2];
02947                      a3[2]=R[cells[i][2]][2]-R[cells[i][6]][2];
02948                      det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
02949                      g1=sqrt(a1[0]*a1[0]+a2[0]*a2[0]+a3[0]*a3[0]);
02950                           g2=sqrt(a1[1]*a1[1]+a2[1]*a2[1]+a3[1]*a3[1]);
02951                           g3=sqrt(a1[2]*a1[2]+a2[2]*a2[2]+a3[2]*a3[2]);
02952                           //need to keep data from previous cell
02953                           if((fabs(det)<eps*eps*eps*det_o)||(det<0)) det=det_o;
02954                           if((fabs(g1)<eps*g1_o)||(g1<0)) g1=g1_o;
02955                           if((fabs(g2)<eps*g2_o)||(g2<0)) g2=g2_o;
02956                           if((fabs(g3)<eps*g3_o)||(g3<0)) g3=g3_o;
02957                           //write to file
02958                           if(me==2) fprintf(stream,"%e %e %e\n0.000000e+00 %e %e\n  0.000000e+00 0.000000e+00 %e\n",1.0/pow(det,1.0/3.0),0.5/pow(det,1.0/3.0),0.5/pow(det,1.0/3.0),0.5*sqrt(3.0)/pow(det,1.0/3.0),0.5/(sqrt(3.0)*pow(det,1.0/3.0)),sqrt(2/3.0)/pow(det,1.0/3.0));
02959               if(me==3) fprintf(stream,"%e %e %e\n0.000000e+00 %e %e\n  0.000000e+00 0.000000e+00 %e\n",1.0/g1,0.5/g2,0.5/g3,0.5*sqrt(3.0)/g2,0.5/(sqrt(3.0)*g3),sqrt(2/3.0)/g3);
02960                           det_o=det; g1_o=g1; g2_o=g2; g3_o=g3;
02961                           Ncells++;
02962                   }
02963           }
02964   }
02965   fclose(stream);
02966 
02967   for(i=0;i<n;i++) free(Q[i]);
02968   free(Q);
02969 
02970   //write new grid connectivity
02971   grid[0]=grid[1]; grid[2]=grid[3];
02972 
02973   stream=fopen(grid,"w+");
02974   fprintf(stream,"%d \n%d \n%d \n0 \n",n,N,Ncells);
02975 
02976   for(i=0;i<N;i++){//node coordinates
02977         for(j=0;j<n;j++) fprintf(stream,"%e ",R[i][j]);
02978         fprintf(stream,"%d \n",mask[i]);
02979   }
02980   for(i=0;i<N;i++) free(R[i]);
02981   free(R);
02982   free(mask);
02983 
02984   for(i=0;i<ncells;i++) if(mcells[i]>=0){//cell connectivity
02985         nvert=0; while(cells[i][nvert]>=0) nvert++;
02986         if((nvert==3)||((n==3)&&(nvert==4))){//tri & tetr
02987                 for(j=0;j<nvert;j++) fprintf(stream,"%d ",cells[i][j]);
02988                 for(j=nvert;j<3*n+n%2;j++) fprintf(stream,"-1 ");
02989             fprintf(stream,"%d \n",mcells[i]);
02990         }
02991         if((n==2)&&(nvert==4)){//quad
02992                 fprintf(stream,"%d %d %d -1 -1 -1 0\n",cells[i][0],cells[i][1],cells[i][2]);
02993         fprintf(stream,"%d %d %d -1 -1 -1 0\n",cells[i][1],cells[i][3],cells[i][0]);
02994         fprintf(stream,"%d %d %d -1 -1 -1 0\n",cells[i][3],cells[i][2],cells[i][1]);
02995         fprintf(stream,"%d %d %d -1 -1 -1 0\n",cells[i][2],cells[i][0],cells[i][3]);
02996         }
02997         if(nvert==8){//hex
02998                 fprintf(stream,"%d %d %d %d -1 -1 -1 -1 -1 -1 0\n",cells[i][0],cells[i][1],cells[i][2],cells[i][4]);
02999         fprintf(stream,"%d %d %d %d -1 -1 -1 -1 -1 -1 0\n",cells[i][1],cells[i][3],cells[i][0],cells[i][5]);
03000         fprintf(stream,"%d %d %d %d -1 -1 -1 -1 -1 -1 0\n",cells[i][3],cells[i][2],cells[i][1],cells[i][7]);
03001         fprintf(stream,"%d %d %d %d -1 -1 -1 -1 -1 -1 0\n",cells[i][2],cells[i][0],cells[i][3],cells[i][6]);
03002                 fprintf(stream,"%d %d %d %d -1 -1 -1 -1 -1 -1 0\n",cells[i][4],cells[i][6],cells[i][5],cells[i][0]);
03003         fprintf(stream,"%d %d %d %d -1 -1 -1 -1 -1 -1 0\n",cells[i][5],cells[i][4],cells[i][7],cells[i][1]);
03004         fprintf(stream,"%d %d %d %d -1 -1 -1 -1 -1 -1 0\n",cells[i][7],cells[i][5],cells[i][6],cells[i][3]);
03005         fprintf(stream,"%d %d %d %d -1 -1 -1 -1 -1 -1 0\n",cells[i][6],cells[i][7],cells[i][4],cells[i][2]);
03006         }
03007   }
03008 
03009 
03010  fclose(stream);
03011 
03012  for(i=0;i<ncells;i++) free(cells[i]);
03013  free(cells);
03014  free(mcells);
03015 
03016 return;
03017 }

double libMesh::VariationalMeshSmoother::minJ	(	int	n,
		int	N,
		LPLPDOUBLE	R,
		LPINT	mask,
		int	ncells,
		LPLPINT	cells,
		LPINT	mcells,
		double	epsilon,
		double	w,
		int	me,
		LPLPLPDOUBLE	H,
		double	vol,
		int	nedges,
		LPINT	edges,
		LPINT	hnodes,
		int	msglev,
		double *	Vmin,
		double *	emax,
		double *	qmin,
		int	adp,
		LPDOUBLE	afun,
		FILE *	sout
	)			`[private]`

Executes one step of minimization algorithm: finds minimization direction (P=H^{-1} J) and solves approximately local minimization problem for optimal step in this minimization direction (tau=min J(R+tau P))

Definition at line 1338 of file mesh_smoother_vsmoother.C.

References std::abs(), alloc_d_n1(), alloc_d_n1_n2(), alloc_d_n1_n2_n3(), alloc_i_n1(), alloc_i_n1_n2(), localP(), std::pow(), and solver().

01342 {
01343 
01344 LPLPLPDOUBLE W; //local Hessian matrix;
01345 LPLPDOUBLE  F, A, G; //F - local gradient; G - adaptation metric;
01346 LPDOUBLE b, u, a; // matrix & rhs for solver, u - solution vector;
01347 LPINT ia, ja; // matrix connectivity for solver;
01348 LPLPINT JA; //A, JA - internal form of global matrix;
01349 LPLPDOUBLE Rpr, P; //P - minimization direction;
01350 double  tau=0.0, J, T, Jpr, lVmin, lemax, lqmin, gVmin=0.0, gemax=0.0,
01351 gqmin=0.0, gtmin0=0.0, gtmax0=0.0, gqmin0=0.0;
01352 double eps, nonzero, Tau_hn, g_i;
01353 int index, i, ii, j, jj, k, l, m, nz;
01354 int sch, ind, columns, nvert, ind_i, ind_j, ind_k;
01355 /* Jpr - value of functional;
01356    nonzero - norm of gradient;
01357    columns - max number of nonzero entries in every row of global matrix;
01358 */
01359 
01360 columns=n*n*10;
01361 W=alloc_d_n1_n2_n3(n,3*n+n%2,3*n+n%2);
01362 F = alloc_d_n1_n2(n,3*n+n%2);
01363 for(i=0;i<n;i++){
01364    F[i]=alloc_d_n1(3*n+n%2);
01365    W[i]=alloc_d_n1_n2(3*n+n%2,3*n+n%2);
01366    for(j=0;j<3*n+n%2;j++) W[i][j]=alloc_d_n1(3*n+n%2);}
01367 Rpr=alloc_d_n1_n2(N,n);
01368 P=alloc_d_n1_n2(N,n);
01369 for(i=0;i<N;i++) {Rpr[i]=alloc_d_n1(n);
01370  P[i]=alloc_d_n1(n);}
01371 A = alloc_d_n1_n2(n*N, columns);
01372 b = alloc_d_n1(n*N);
01373 u = alloc_d_n1(n*N);
01374 a = alloc_d_n1(n*N*columns);
01375 ia = alloc_i_n1(n*N+1);
01376 ja = alloc_i_n1(n*N*columns);
01377 JA = alloc_i_n1_n2(n*N, columns);
01378 for(i=0;i<n*N;i++) {
01379    A[i] = alloc_d_n1(columns);
01380    JA[i] = alloc_i_n1(columns);}
01381 G=alloc_d_n1_n2(ncells,n);
01382 for(i=0;i<ncells;i++) G[i]=alloc_d_n1(n);
01383 
01384 //---------find minimization direction P-----------------
01385 nonzero=0; Jpr=0;
01386 for(i=0; i<n*N; i++){//initialise matrix and rhs
01387     for(j=0; j<columns; j++){ A[i][j] = 0; JA[i][j] =0;}
01388     b[i] = 0;
01389     }
01390 for(i=0; i<ncells; i++)
01391 {
01392   nvert=0;
01393   while(cells[i][nvert]>=0)
01394     nvert++;
01395   //---determination of local matrices on each cell----
01396 
01397   for(j=0;j<n;j++)
01398   {
01399     G[i][j]=0; //adaptation metric G is held constant throughout minJ run
01400     if(adp<0)
01401     {
01402       for(k=0;k<abs(adp);k++)
01403         G[i][j]=G[i][j]+afun[i*(-adp)+k]; //cell-based adaptivity is computed here
01404     }
01405   }
01406      for(index=0;index<n;index++){//initialise local matrices
01407          for(k=0;k<3*n+n%2;k++){ F[index][k]=0;
01408              for(j=0;j<3*n+n%2;j++) W[index][k][j]=0;
01409          }
01410      }
01411      if(mcells[i]>=0){//if cell is not excluded
01412          Jpr+=localP(n, W, F, R, cells[i], mask, epsilon, w, nvert, H[i],
01413                      me, vol, 0, &lVmin, &lqmin, adp, afun, G[i], sout);
01414      } else {
01415             for(index=0;index<n;index++) for(j=0;j<nvert;j++) W[index][j][j]=1;
01416      }
01417     //----assembly of an upper triangular part of a global matrix A------
01418     for(index=0;index<n;index++)
01419     {
01420       for(l=0; l<nvert; l++)
01421       {
01422         for(m=0; m<nvert; m++){
01423           if((W[index][l][m]!=0)&&(cells[i][m]>=cells[i][l]))
01424           {
01425             sch=0; ind=1;
01426             while (ind!=0)
01427             {
01428               if(A[cells[i][l]+index*N][sch]!=0)
01429               {
01430                   if(JA[cells[i][l]+index*N][sch]==(cells[i][m]+index*N))
01431                   {
01432                     A[cells[i][l]+index*N][sch] = A[cells[i][l]+index*N][sch] + W[index][l][m];
01433                     ind=0;
01434                   }
01435                   else
01436                     sch++;
01437               }
01438               else
01439               {
01440                 A[cells[i][l]+index*N][sch] = W[index][l][m];
01441                 JA[cells[i][l]+index*N][sch]=cells[i][m]+index*N;
01442                 ind = 0;
01443               }
01444               if(sch>columns-1) fprintf(sout,"error: # of nonzero entries in the %d row of Hessian =%d >= columns=%d \n",cells[i][l],sch,columns);
01445             }
01446           }
01447         }
01448         b[cells[i][l]+index*N]=b[cells[i][l]+index*N]-F[index][l];
01449       }
01450     }
01451 //------end of matrix A---------
01452   }
01453  //-------HN correction--------
01454  Tau_hn=pow(vol,1.0/(double)n)*1e-10; //tolerance for HN being the mid-edge node for its parents
01455  for(i=0;i<nedges;i++){
01456      ind_i=hnodes[i]; ind_j=edges[2*i]; ind_k=edges[2*i+1];
01457      for(j=0;j<n;j++){
01458          g_i=R[ind_i][j]-0.5*(R[ind_j][j]+R[ind_k][j]);
01459          Jpr+=g_i*g_i/(2*Tau_hn);
01460          b[ind_i+j*N]-=g_i/Tau_hn;
01461          b[ind_j+j*N]+=0.5*g_i/Tau_hn;
01462          b[ind_k+j*N]+=0.5*g_i/Tau_hn;
01463      }
01464      for(ind=0;ind<columns;ind++){
01465          if(JA[ind_i][ind]==ind_i) A[ind_i][ind]+=1.0/Tau_hn;
01466          if(JA[ind_i+N][ind]==ind_i+N) A[ind_i+N][ind]+=1.0/Tau_hn;
01467          if(n==3) if(JA[ind_i+2*N][ind]==ind_i+2*N) A[ind_i+2*N][ind]+=1.0/Tau_hn;
01468          if((JA[ind_i][ind]==ind_j)||(JA[ind_i][ind]==ind_k)) A[ind_i][ind]-=0.5/Tau_hn;
01469          if((JA[ind_i+N][ind]==ind_j+N)||(JA[ind_i+N][ind]==ind_k+N)) A[ind_i+N][ind]-=0.5/Tau_hn;
01470          if(n==3) if((JA[ind_i+2*N][ind]==ind_j+2*N)||(JA[ind_i+2*N][ind]==ind_k+2*N)) A[ind_i+2*N][ind]-=0.5/Tau_hn;
01471          if(JA[ind_j][ind]==ind_i) A[ind_j][ind]-=0.5/Tau_hn;
01472          if(JA[ind_k][ind]==ind_i) A[ind_k][ind]-=0.5/Tau_hn;
01473          if(JA[ind_j+N][ind]==ind_i+N) A[ind_j+N][ind]-=0.5/Tau_hn;
01474          if(JA[ind_k+N][ind]==ind_i+N) A[ind_k+N][ind]-=0.5/Tau_hn;
01475          if(n==3) if(JA[ind_j+2*N][ind]==ind_i+2*N) A[ind_j+2*N][ind]-=0.5/Tau_hn;
01476          if(n==3) if(JA[ind_k+2*N][ind]==ind_i+2*N) A[ind_k+2*N][ind]-=0.5/Tau_hn;
01477          if((JA[ind_j][ind]==ind_j)||(JA[ind_j][ind]==ind_k)) A[ind_j][ind]+=0.25/Tau_hn;
01478          if((JA[ind_k][ind]==ind_j)||(JA[ind_k][ind]==ind_k)) A[ind_k][ind]+=0.25/Tau_hn;
01479          if((JA[ind_j+N][ind]==ind_j+N)||(JA[ind_j+N][ind]==ind_k+N)) A[ind_j+N][ind]+=0.25/Tau_hn;
01480          if((JA[ind_k+N][ind]==ind_j+N)||(JA[ind_k+N][ind]==ind_k+N)) A[ind_k+N][ind]+=0.25/Tau_hn;
01481          if(n==3) if((JA[ind_j+2*N][ind]==ind_j+2*N)||(JA[ind_j+2*N][ind]==ind_k+2*N)) A[ind_j+2*N][ind]+=0.25/Tau_hn;
01482          if(n==3) if((JA[ind_k+2*N][ind]==ind_j+2*N)||(JA[ind_k+2*N][ind]==ind_k+2*N)) A[ind_k+2*N][ind]+=0.25/Tau_hn;
01483      }
01484  }
01485 //------||\grad J||_2------
01486 for(i=0;i<n*N;i++){ nonzero=nonzero+b[i]*b[i];}
01487 
01488 //-----sort matrix A-----
01489   for(i=0;i<n*N;i++)
01490   {
01491       for(j=columns-1;j>1;j--)
01492       {
01493           for(k=0;k<j;k++)
01494           {
01495              if(JA[i][k]>JA[i][k+1])
01496              {
01497                sch=JA[i][k];
01498                JA[i][k]=JA[i][k+1];
01499                JA[i][k+1]=sch;
01500                tau=A[i][k];
01501                A[i][k]=A[i][k+1];
01502                A[i][k+1]=tau;
01503              }
01504          }
01505       }
01506   }
01507 
01508 eps=sqrt(vol)*1e-9;
01509 
01510 //-------solver for P (unconstrained)-------
01511   ia[0]=0; j=0;
01512   for(i=0;i<n*N;i++)
01513   {
01514       u[i]=0;
01515       nz=0;
01516       for(k=0;k<columns;k++){
01517          if(A[i][k]!=0)
01518          {
01519            nz++;
01520            ja[j]=JA[i][k]+1;
01521            a[j]=A[i][k];
01522            j++;
01523          }
01524       }
01525       ia[i+1]=ia[i]+nz;
01526   }
01527   nz=ia[n*N];
01528   m=n*N;
01529   if(msglev>=3) sch=1; else sch=0;
01530 
01531 
01532   nz=solver(m,ia,ja,a,u,b,eps,100,sch, sout);
01533 //    sol_pcg_pj(m,ia,ja,a,u,b,eps,100,sch);
01534 
01535   for(i=0;i<N;i++){ //ensure fixed nodes are not moved
01536       for(index=0;index<n;index++)
01537           if(mask[i]==1) P[i][index]=0; else P[i][index]=u[i+index*N];
01538   }
01539  //----P is determined--------
01540  if(msglev>=4){
01541  for(i=0;i<N;i++){
01542     for(j=0;j<n;j++)  if(P[i][j]!=0) fprintf(sout, "P[%d][%d]=%f  ",i,j,P[i][j]);
01543     }
01544  }
01545 //-------local minimization problem, determination of tau----------
01546 if(msglev>=3) fprintf(sout, "dJ=%e J0=%e \n",sqrt(nonzero),Jpr);
01547 J=1e32; j=1;
01548 while((Jpr<=J)&&(j>-30)){
01549    j=j-1;
01550    tau=pow(2.0,j); //search through finite # of values for tau
01551    for(i=0;i<N;i++)
01552    {
01553        for(k=0;k<n;k++)
01554          Rpr[i][k]=R[i][k]+tau*P[i][k];
01555    }
01556    J=0;
01557    gVmin=1e32; gemax=-1e32; gqmin=1e32;
01558    for(i=0; i<ncells; i++)
01559    {
01560      if(mcells[i]>=0)
01561      {
01562         nvert=0;
01563         while(cells[i][nvert]>=0)
01564           nvert++;
01565         lemax=localP(n, W, F, Rpr, cells[i], mask, epsilon, w, nvert, H[i], me, vol, 1, &lVmin, &lqmin, adp, afun, G[i], sout);
01566         J+=lemax;
01567         if(gVmin>lVmin)
01568           gVmin=lVmin;
01569         if(gemax<lemax)
01570           gemax=lemax;
01571         if(gqmin>lqmin)
01572           gqmin=lqmin;
01573 
01574         //------HN correction--------
01575         for(ii=0;ii<nedges;ii++)
01576         {
01577             ind_i=hnodes[ii];
01578             ind_j=edges[2*ii];
01579             ind_k=edges[2*ii+1];
01580             for(jj=0;jj<n;jj++)
01581             {
01582                 g_i=Rpr[ind_i][jj]-0.5*(Rpr[ind_j][jj]+Rpr[ind_k][jj]);
01583                 J+=g_i*g_i/(2*Tau_hn);
01584             }
01585          }
01586       }
01587    }
01588    if(msglev>=3)
01589      fprintf(sout, "tau=%f J=%f \n",tau,J);
01590 }
01591 
01592 if(j==-30)
01593   T=0;
01594 else
01595 {
01596   Jpr=J;
01597   for(i=0;i<N;i++)
01598   {
01599     for(k=0;k<n;k++)
01600       Rpr[i][k]=R[i][k]+tau*0.5*P[i][k];
01601   }
01602   J=0; gtmin0=1e32; gtmax0=-1e32; gqmin0=1e32;
01603   for(i=0; i<ncells; i++)
01604   {
01605     if(mcells[i]>=0)
01606     {
01607       nvert=0; while(cells[i][nvert]>=0) nvert++;
01608       lemax=localP(n, W, F, Rpr, cells[i], mask, epsilon, w, nvert, H[i], me, vol, 1, &lVmin, &lqmin, adp, afun, G[i], sout);
01609       J+=lemax;
01610       if(gtmin0>lVmin)
01611         gtmin0=lVmin;
01612       if(gtmax0<lemax)
01613         gtmax0=lemax;
01614       if(gqmin0>lqmin)
01615         gqmin0=lqmin;
01616       //-------HN correction--------
01617        for(ii=0;ii<nedges;ii++){
01618            ind_i=hnodes[ii]; ind_j=edges[2*ii]; ind_k=edges[2*ii+1];
01619            for(jj=0;jj<n;jj++){
01620                g_i=Rpr[ind_i][jj]-0.5*(Rpr[ind_j][jj]+Rpr[ind_k][jj]);
01621                J+=g_i*g_i/(2*Tau_hn);
01622            }
01623        }//HN
01624    }
01625   }
01626 }
01627  if(Jpr>J) {
01628     T=0.5*tau;
01629     (*Vmin)=gtmin0;
01630     (*emax)=gtmax0;
01631     (*qmin)=gqmin0;
01632  } else {
01633     T=tau; J=Jpr;
01634     (*Vmin)=gVmin;
01635     (*emax)=gemax;
01636     (*qmin)=gqmin;
01637  }
01638 
01639  nonzero=0;
01640 for(j=0;j<N;j++) for(k=0;k<n;k++){
01641      R[j][k]=R[j][k]+T*P[j][k];
01642      nonzero+=T*P[j][k]*T*P[j][k];
01643 }
01644 
01645 if(msglev>=2)
01646   fprintf(sout, "tau=%e, J=%e  \n",T,J);
01647 
01648 free(b);
01649 free(u);
01650 free(ia);
01651 free(ja);
01652 free(a);
01653 for(i=0;i<n;i++){
01654    for(j=0;j<3*n+n%2;j++) free(W[i][j]);
01655    free(W[i]);
01656    free(F[i]);}
01657 free(W);
01658 free(F);
01659 for(i=0;i<N;i++) {free(Rpr[i]);
01660     free(P[i]);}
01661 free(Rpr);
01662 free(P);
01663 for(i=0;i<n*N;i++) {
01664     free(A[i]);
01665     free(JA[i]);}
01666 free(A);
01667 free(JA);
01668 for(i=0;i<ncells;i++) free(G[i]);
01669 free(G);
01670 
01671 return (sqrt(nonzero));
01672 
01673 }

double libMesh::VariationalMeshSmoother::minJ_BC	(	int	N,
		LPLPDOUBLE	R,
		LPINT	mask,
		int	ncells,
		LPLPINT	cells,
		LPINT	mcells,
		double	epsilon,
		double	w,
		int	me,
		LPLPLPDOUBLE	H,
		double	vol,
		int	msglev,
		double *	Vmin,
		double *	emax,
		double *	qmin,
		int	adp,
		LPDOUBLE	afun,
		int	NCN,
		FILE *	sout
	)			`[private]`

minJ() with sliding Boundary Nodes constraints and no account for HN, using Lagrange multiplier formulation: minimize L=J+ lam*g; only works in 2D

Definition at line 1680 of file mesh_smoother_vsmoother.C.

References std::abs(), alloc_d_n1(), alloc_d_n1_n2(), alloc_d_n1_n2_n3(), alloc_i_n1(), localP(), and std::pow().

01683 {
01684 /*---new form of matrices, 5 iterations for minL---*/
01685 LPLPLPDOUBLE W;
01686 LPLPDOUBLE  F, G;
01687 LPDOUBLE b, hm, Plam, constr, lam;
01688 LPINT Bind;
01689 LPLPDOUBLE Rpr, P;
01690 double  tau=0.0, J=0.0, T, Jpr, L, lVmin, lqmin, gVmin=0.0, gqmin=0.0, gVmin0=0.0,
01691 gqmin0=0.0, lemax, gemax=0.0, gemax0=0.0;
01692 double a, g, qq=0.0, eps, nonzero, x0, y0, del1, del2, Bx, By;
01693 int index, i, j, k=0, l, nz, I;
01694 int ind, nvert;
01695 
01696 //----memory------
01697 Bind=alloc_i_n1(NCN); //array of sliding BN
01698 lam=alloc_d_n1(NCN);
01699 hm=alloc_d_n1(2*N);
01700 Plam=alloc_d_n1(NCN);
01701 constr=alloc_d_n1(4*NCN); //holds constraints = local approximation to the boundary
01702 F = alloc_d_n1_n2(2, 6);
01703 W=alloc_d_n1_n2_n3(2, 6, 6);
01704 for(i=0;i<2;i++){
01705    F[i]=alloc_d_n1(6);
01706    W[i]=alloc_d_n1_n2(6,6);
01707    for(j=0;j<6;j++) W[i][j]=alloc_d_n1(6);}
01708 Rpr=alloc_d_n1_n2(N,2);
01709 P=alloc_d_n1_n2(N,2);
01710 for(i=0;i<N;i++) {Rpr[i]=alloc_d_n1(2);
01711  P[i]=alloc_d_n1(2);}
01712 b = alloc_d_n1(2*N);
01713 G=alloc_d_n1_n2(ncells,6);
01714 for(i=0;i<ncells;i++) G[i]=alloc_d_n1(6);
01715 
01716 //-----assembler of constraints-----
01717   eps=sqrt(vol)*1e-9;
01718   for(i=0;i<4*NCN;i++) constr[i]=1.0/eps;
01719   I=0; for(i=0;i<N;i++) if(mask[i]==2) {Bind[I]=i; I++;}
01720   for(I=0;I<NCN;I++){ i=Bind[I];
01721       ind=0;
01722       //---boundary connectivity----
01723       for(l=0;l<ncells;l++){
01724           nvert=0; while(cells[l][nvert]>=0) nvert++;
01725                   switch(nvert)
01726                   {case 3: //tri
01727                   if(i==cells[l][0]){
01728                           if(mask[cells[l][1]]>0) {if(ind==0) {j=cells[l][1]; ind++;} else {k=cells[l][1]; ind++;}}
01729                           if(mask[cells[l][2]]>0) {if(ind==0) {j=cells[l][2]; ind++;} else {k=cells[l][2]; ind++;}}
01730                   }
01731                   if(i==cells[l][1]){
01732                           if(mask[cells[l][0]]>0) {if(ind==0) {j=cells[l][0]; ind++;} else {k=cells[l][0]; ind++;}}
01733                           if(mask[cells[l][2]]>0) {if(ind==0) {j=cells[l][2]; ind++;} else {k=cells[l][2]; ind++;}}
01734                   }
01735                   if(i==cells[l][2]){
01736                           if(mask[cells[l][1]]>0) {if(ind==0) {j=cells[l][1]; ind++;} else {k=cells[l][1]; ind++;}}
01737                           if(mask[cells[l][0]]>0) {if(ind==0) {j=cells[l][0]; ind++;} else {k=cells[l][0]; ind++;}}
01738                   }
01739                   break;
01740                   case 4: //quad
01741                   if((i==cells[l][0])||(i==cells[l][3])){
01742                           if(mask[cells[l][1]]>0) {if(ind==0) {j=cells[l][1]; ind++;} else {k=cells[l][1]; ind++;}}
01743                           if(mask[cells[l][2]]>0) {if(ind==0) {j=cells[l][2]; ind++;} else {k=cells[l][2]; ind++;}}
01744                   }
01745                   if((i==cells[l][1])||(i==cells[l][2])){
01746                           if(mask[cells[l][0]]>0) {if(ind==0) {j=cells[l][0]; ind++;} else {k=cells[l][0]; ind++;}}
01747                           if(mask[cells[l][3]]>0) {if(ind==0) {j=cells[l][3]; ind++;} else {k=cells[l][3]; ind++;}}
01748                   }
01749                   break;
01750                   case 6: //quad tri
01751                   if(i==cells[l][0]){
01752                           if(mask[cells[l][1]]>0) {if(ind==0) {j=cells[l][5]; ind++;} else {k=cells[l][5]; ind++;}}
01753                           if(mask[cells[l][2]]>0) {if(ind==0) {j=cells[l][4]; ind++;} else {k=cells[l][4]; ind++;}}
01754                   }
01755                   if(i==cells[l][1]){
01756                           if(mask[cells[l][0]]>0) {if(ind==0) {j=cells[l][5]; ind++;} else {k=cells[l][5]; ind++;}}
01757                           if(mask[cells[l][2]]>0) {if(ind==0) {j=cells[l][3]; ind++;} else {k=cells[l][3]; ind++;}}
01758                   }
01759                   if(i==cells[l][2]){
01760                           if(mask[cells[l][1]]>0) {if(ind==0) {j=cells[l][3]; ind++;} else {k=cells[l][3]; ind++;}}
01761                           if(mask[cells[l][0]]>0) {if(ind==0) {j=cells[l][4]; ind++;} else {k=cells[l][4]; ind++;}}
01762                   }
01763                   if(i==cells[l][3]){j=cells[l][1]; k=cells[l][2];}
01764           if(i==cells[l][4]){j=cells[l][2]; k=cells[l][0];}
01765           if(i==cells[l][5]){j=cells[l][0]; k=cells[l][1];}
01766                   break;}
01767           }//end boundary connectivity
01768       //lines
01769       del1=R[j][0]-R[i][0]; del2=R[i][0]-R[k][0];
01770       if((fabs(del1)<eps)&&(fabs(del2)<eps)) {constr[I*4]=1; constr[I*4+1]=0;
01771           constr[I*4+2]=R[i][0]; constr[I*4+3]=R[i][1];} else {
01772           del1=R[j][1]-R[i][1]; del2=R[i][1]-R[k][1];
01773           if((fabs(del1)<eps)&&(fabs(del2)<eps)) {constr[I*4]=0; constr[I*4+1]=1;
01774               constr[I*4+2]=R[i][0]; constr[I*4+3]=R[i][1];
01775               } else {
01776                   J=(R[i][0]-R[j][0])*(R[k][1]-R[j][1])-(R[k][0]-R[j][0])*(R[i][1]-R[j][1]);
01777                   if(fabs(J)<eps) {
01778                       constr[I*4]=1.0/(R[k][0]-R[j][0]); constr[I*4+1]=-1.0/(R[k][1]-R[j][1]);
01779                       constr[I*4+2]=R[i][0]; constr[I*4+3]=R[i][1];
01780                   } else { //circle
01781                       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])+
01782                           (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);
01783                       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])+
01784                           (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);
01785                       constr[I*4]=x0; constr[I*4+1]=y0;
01786                       constr[I*4+2]=(R[i][0]-x0)*(R[i][0]-x0)+(R[i][1]-y0)*(R[i][1]-y0);
01787                   }
01788               }
01789       }
01790   }
01791 /*for(i=0;i<NCN;i++){
01792     for(j=0;j<4;j++) fprintf(stdout," %e ",constr[4*i+j]);
01793     fprintf(stdout, "constr %d node %d \n",i,Bind[i]);
01794 }*/
01795 
01796 
01797 
01798 //----initial values----
01799 for(i=0;i<NCN;i++) lam[i]=0;
01800 for(nz=0;nz<5;nz++){
01801   //---------find H and -grad J-----------------
01802    nonzero=0; Jpr=0;
01803    for(i=0; i<2*N; i++){ b[i] = 0; hm[i]=0; }
01804    for(i=0; i<ncells; i++){
01805         nvert=0; while(cells[i][nvert]>=0) nvert++;
01806         for(j=0;j<nvert;j++){ G[i][j]=0;
01807         if(adp<0) for(k=0;k<abs(adp);k++) G[i][j]=G[i][j]+afun[i*(-adp)+k]; }
01808         for(index=0;index<2;index++){
01809             for(k=0;k<nvert;k++){ F[index][k]=0;
01810                 for(j=0;j<nvert;j++) W[index][k][j]=0;
01811             }
01812         }
01813         if(mcells[i]>=0){
01814                      Jpr+=localP(2, W, F, R, cells[i], mask, epsilon, w, nvert, H[i],
01815                          me, vol, 0, &lVmin, &lqmin, adp, afun, G[i], sout);
01816         } else {
01817             for(index=0;index<2;index++) for(j=0;j<nvert;j++) W[index][j][j]=1;
01818         }
01819         for(index=0;index<2;index++){
01820             for(l=0; l<nvert; l++){//-diagonal Hessian-
01821                 hm[cells[i][l]+index*N]=hm[cells[i][l]+index*N]+W[index][l][l];
01822                 b[cells[i][l]+index*N]=b[cells[i][l]+index*N]-F[index][l];
01823             }
01824         }
01825     }
01826     //------||grad J||_2------
01827     for(i=0;i<2*N;i++){ nonzero=nonzero+b[i]*b[i];}
01828     //-----solve for Plam--------
01829     for(I=0;I<NCN;I++){i=Bind[I];
01830        if(constr[4*I+3]<0.5/eps) {Bx=constr[4*I]; By=constr[4*I+1];
01831            g=(R[i][0]-constr[4*I+2])*constr[4*I]+(R[i][1]-constr[4*I+3])*constr[4*I+1];
01832        } else {
01833           Bx=2*(R[i][0]-constr[4*I]); By=2*(R[i][1]-constr[4*I+1]);
01834               hm[i]+=2*lam[I]; hm[i+N]+=2*lam[I];
01835               g=(R[i][0]-constr[4*I])*(R[i][0]-constr[4*I])+(R[i][1]-constr[4*I+1])*(R[i][1]-constr[4*I+1])-constr[4*I+2];
01836        }
01837        Jpr+=lam[I]*g;
01838        qq=Bx*b[i]/hm[i]+By*b[i+N]/hm[i+N]-g;
01839        a=Bx*Bx/hm[i]+By*By/hm[i+N];
01840        if(a!=0) Plam[I]=qq/a; else fprintf(sout,"error: B^TH-1B is degenerate \n");
01841        b[i]-=Plam[I]*Bx; b[i+N]-=Plam[I]*By;
01842        Plam[I]-=lam[I];
01843    }
01844    //-----------solve for P------------
01845    for(i=0;i<N;i++) {P[i][0]=b[i]/hm[i]; P[i][1]=b[i+N]/hm[i+N];}
01846    //-----correct solution-----
01847    for(i=0;i<N;i++) for(j=0;j<2;j++) if((fabs(P[i][j])<eps)||(mask[i]==1)) P[i][j]=0;
01848    //----P is determined--------
01849    if(msglev>=3){
01850        for(i=0;i<N;i++){
01851            for(j=0;j<2;j++)  if(P[i][j]!=0) fprintf(sout, "P[%d][%d]=%f  ",i,j,P[i][j]);
01852        }
01853    }
01854    //-------local minimization problem, determination of tau----------
01855    if(msglev>=3) fprintf(sout, "dJ=%e L0=%e \n",sqrt(nonzero), Jpr);
01856    L=1e32; j=1;
01857    while((Jpr<=L)&&(j>-30)){
01858       j=j-1;
01859       tau=pow(2.0,j);
01860       for(i=0;i<N;i++){
01861           for(k=0;k<2;k++) Rpr[i][k]=R[i][k]+tau*P[i][k];}
01862       J=0; gVmin=1e32; gemax=-1e32; gqmin=1e32;
01863       for(i=0; i<ncells; i++) if(mcells[i]>=0){
01864           nvert=0; while(cells[i][nvert]>=0) nvert++;
01865                   lemax=localP(2, W, F, Rpr, cells[i], mask, epsilon, w, nvert, H[i], me, vol, 1, &lVmin,
01866                        &lqmin, adp, afun, G[i], sout);
01867               J+=lemax;
01868           if(gVmin>lVmin) gVmin=lVmin;
01869           if(gemax<lemax) gemax=lemax;
01870           if(gqmin>lqmin) gqmin=lqmin;
01871           }
01872       L=J;
01873       //----constraints contribution----
01874       for(I=0;I<NCN;I++){i=Bind[I];
01875           if(constr[4*I+3]<0.5/eps) g=(Rpr[i][0]-constr[4*I+2])*constr[4*I]+(Rpr[i][1]-constr[4*I+3])*constr[4*I+1];
01876               else g=(Rpr[i][0]-constr[4*I])*(Rpr[i][0]-constr[4*I])+
01877                    (Rpr[i][1]-constr[4*I+1])*(Rpr[i][1]-constr[4*I+1])-constr[4*I+2];
01878           L+=(lam[I]+tau*Plam[I])*g;
01879           }
01880       //----end of constraints----
01881       if(msglev>=3) fprintf(sout," tau=%f J=%f \n",tau,J);
01882    }
01883    if(j==-30) { T=0; } else {
01884       Jpr=L; qq=J;
01885       for(i=0;i<N;i++){
01886           for(k=0;k<2;k++) Rpr[i][k]=R[i][k]+tau*0.5*P[i][k];}
01887       J=0; gVmin0=1e32; gemax0=-1e32; gqmin0=1e32;
01888       for(i=0; i<ncells; i++) if(mcells[i]>=0){
01889           nvert=0; while(cells[i][nvert]>=0) nvert++;
01890                   lemax=localP(2, W, F, Rpr, cells[i], mask, epsilon, w, nvert, H[i], me, vol, 1, &lVmin,
01891                        &lqmin, adp, afun, G[i], sout);
01892                   J+=lemax;
01893               if(gVmin0>lVmin) gVmin0=lVmin;
01894           if(gemax0<lemax) gemax0=lemax;
01895           if(gqmin0>lqmin) gqmin0=lqmin;
01896           }
01897       L=J;
01898       //----constraints contribution----
01899       for(I=0;I<NCN;I++){i=Bind[I];
01900           if(constr[4*I+3]<0.5/eps) g=(Rpr[i][0]-constr[4*I+2])*constr[4*I]+(Rpr[i][1]-constr[4*I+3])*constr[4*I+1];
01901               else g=(Rpr[i][0]-constr[4*I])*(Rpr[i][0]-constr[4*I])+
01902                      (Rpr[i][1]-constr[4*I+1])*(Rpr[i][1]-constr[4*I+1])-constr[4*I+2];
01903           L+=(lam[I]+tau*0.5*Plam[I])*g;
01904           }
01905       //----end of constraints----
01906    }
01907    if(Jpr>L) {
01908       T=0.5*tau;
01909       (*Vmin)=gVmin0;
01910       (*emax)=gemax0;
01911       (*qmin)=gqmin0;
01912    } else {
01913       T=tau; L=Jpr; J=qq;
01914       (*Vmin)=gVmin;
01915       (*emax)=gemax;
01916       (*qmin)=gqmin;
01917    }
01918 
01919    for(j=0;j<N;j++) for(k=0;k<2;k++) R[j][k]+=T*P[j][k];
01920    for(i=0;i<NCN;i++) lam[i]+=T*Plam[i];
01921 
01922 }//end Lagrangian iter
01923 
01924  if(msglev>=2) fprintf(sout, "tau=%e, J=%e, L=%e  \n",T,J,L);
01925 
01926 free(lam);
01927 free(b);
01928 for(i=0;i<2;i++){
01929    for(j=0;j<6;j++) free(W[i][j]);
01930    free(W[i]);
01931    free(F[i]);}
01932 free(W);
01933 free(F);
01934 for(i=0;i<N;i++) {free(Rpr[i]);
01935     free(P[i]);}
01936 free(Rpr);
01937 free(P);
01938 for(i=0;i<ncells;i++) free(G[i]);
01939 free(G);
01940 free(Bind);
01941 free(constr);
01942 free(hm);
01943 free(Plam);
01944 
01945 return (sqrt(nonzero));
01946 
01947 }

double libMesh::VariationalMeshSmoother::minJ_l	(	int	n,
		int	N,
		LPLPDOUBLE	R,
		LPINT	mask,
		int	ncells,
		LPLPINT	cells,
		LPINT	mcells,
		double	epsilon,
		double	w,
		int	me,
		LPLPLPDOUBLE	H,
		double	vol,
		int	nedges,
		LPINT	edges,
		LPINT	hnodes,
		int	msglev,
		double *	Vmin,
		double *	emax,
		double *	qmin,
		int	adp,
		LPDOUBLE	afun,
		FILE *	sout
	)			`[private]`

double libMesh::VariationalMeshSmoother::minq	(	int	n,
		int	N,
		LPLPDOUBLE	R,
		int	ncells,
		LPLPINT	cells,
		LPINT	mcells,
		int	me,
		LPLPLPDOUBLE	H,
		double *	vol,
		double *	Vmin,
		FILE *	sout
	)			`[private]`

Compute min Jacobian determinant (minq), min cell volume (Vmin), and average cell volume (vol).

Definition at line 1154 of file mesh_smoother_vsmoother.C.

References alloc_d_n1(), alloc_d_n1_n2(), basisA(), jac2(), and jac3().

01156 {
01157   LPLPDOUBLE Q;
01158   double K[9], a1[3], a2[3], a3[3];
01159   double  v, vmin, gqmin, det, vcell, sigma=0.0;
01160   int  ii, i, j, k, l, m, nn, kk, ll;
01161 
01162   Q = alloc_d_n1_n2(3, 10);
01163   for(i=0;i<n;i++) Q[i]=alloc_d_n1(3*n+n%2);
01164 
01165   v=0; vmin=1e32; gqmin=1e32;
01166 
01167 for(ii=0; ii<ncells; ii++) if(mcells[ii]>=0){
01168     if(n==2){//2D
01169         if(cells[ii][3]==-1){//tri
01170           basisA(2,Q,3,K,H[ii],me);
01171           for(k=0;k<2;k++){a1[k]=0; a2[k]=0;
01172              for(l=0;l<3;l++){
01173                  a1[k]=a1[k]+Q[k][l]*R[cells[ii][l]][0];
01174                  a2[k]=a2[k]+Q[k][l]*R[cells[ii][l]][1];
01175              }
01176           }
01177           det=jac2(a1[0],a1[1],a2[0],a2[1]);
01178           if(gqmin>det) gqmin=det;
01179           if(vmin>det) vmin=det;
01180           v+=det;
01181         } else if(cells[ii][4]==-1){ vcell=0; //quad
01182           for(i=0; i<2; i++){ K[0]=i;
01183             for(j=0; j<2; j++){ K[1]=j;
01184                basisA(2,Q,4,K,H[ii],me);
01185                for(k=0;k<2;k++){a1[k]=0; a2[k]=0;
01186                   for(l=0;l<4;l++){
01187                       a1[k]=a1[k]+Q[k][l]*R[cells[ii][l]][0];
01188                       a2[k]=a2[k]+Q[k][l]*R[cells[ii][l]][1];
01189                   }
01190                }
01191                det=jac2(a1[0],a1[1],a2[0],a2[1]);
01192                if(gqmin>det) gqmin=det;
01193                            v+=0.25*det;
01194                            vcell+=0.25*det;
01195             }
01196           }
01197           if(vmin>vcell) vmin=vcell;
01198                 } else { vcell=0; //quad tri
01199             for(i=0; i<3; i++){ K[0]=i*0.5; k=i/2; K[1]=(double)k;
01200                for(j=0; j<3; j++){  K[2]=j*0.5; k=j/2; K[3]=(double)k;
01201                               basisA(2,Q,6,K,H[ii],me);
01202                   for(k=0;k<2;k++){a1[k]=0; a2[k]=0;
01203                      for(l=0;l<6;l++){
01204                         a1[k]=a1[k]+Q[k][l]*R[cells[ii][l]][0];
01205                         a2[k]=a2[k]+Q[k][l]*R[cells[ii][l]][1];
01206                      }
01207                   }
01208                   det=jac2(a1[0],a1[1],a2[0],a2[1]);
01209                   if(gqmin>det) gqmin=det;
01210                                   if(i==j) sigma=1.0/12; else sigma=1.0/24;
01211                                   v+=sigma*det;
01212                                   vcell+=sigma*det;
01213                            }
01214                         }
01215                         if(vmin>vcell) vmin=vcell;
01216                 }
01217         }
01218     if(n==3){//3D
01219        if(cells[ii][4]==-1){//tetr
01220            basisA(3,Q,4,K,H[ii],me);
01221            for(k=0;k<3;k++){a1[k]=0; a2[k]=0; a3[k]=0;
01222               for(l=0;l<4;l++){
01223                   a1[k]=a1[k]+Q[k][l]*R[cells[ii][l]][0];
01224                   a2[k]=a2[k]+Q[k][l]*R[cells[ii][l]][1];
01225                   a3[k]=a3[k]+Q[k][l]*R[cells[ii][l]][2];
01226               }
01227            }
01228            det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
01229            if(gqmin>det) gqmin=det;
01230                    if(vmin>det) vmin=det;
01231            v+=det;
01232        } else if(cells[ii][6]==-1){//prism
01233           vcell=0;
01234           for(i=0;i<2;i++){ K[0]=i;
01235               for(j=0;j<2;j++){ K[1]=j;
01236                   for(k=0;k<3;k++) {K[2]=(double)k/2.0; K[3]=(double)(k%2);
01237                       basisA(3,Q,6,K,H[ii],me);
01238                       for(kk=0;kk<3;kk++){a1[kk]=0; a2[kk]=0; a3[kk]=0;
01239                           for(ll=0;ll<6;ll++){
01240                               a1[kk]=a1[kk]+Q[kk][ll]*R[cells[ii][ll]][0];
01241                               a2[kk]=a2[kk]+Q[kk][ll]*R[cells[ii][ll]][1];
01242                               a3[kk]=a3[kk]+Q[kk][ll]*R[cells[ii][ll]][2];
01243                           }
01244                       }
01245                       det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
01246                       if(gqmin>det) gqmin=det; sigma=1.0/12.0;
01247                       v+=sigma*det; vcell+=sigma*det;
01248                                   }
01249                           }
01250                   }
01251           if(vmin>vcell) vmin=vcell;
01252            } else if(cells[ii][8]==-1){ vcell=0; //hex
01253          for(i=0; i<2; i++){ K[0]=i;
01254              for(j=0; j<2; j++){ K[1]=j;
01255                  for(k=0; k<2; k++){ K[2]=k;
01256                      for(l=0; l<2; l++){ K[3]=l;
01257                          for(m=0; m<2; m++){ K[4]=m;
01258                              for(nn=0; nn<2; nn++){ K[5]=nn;
01259                                 basisA(3,Q,8,K,H[ii],me);
01260                                 for(kk=0;kk<3;kk++){a1[kk]=0; a2[kk]=0; a3[kk]=0;
01261                                    for(ll=0;ll<8;ll++){
01262                                       a1[kk]=a1[kk]+Q[kk][ll]*R[cells[ii][ll]][0];
01263                                       a2[kk]=a2[kk]+Q[kk][ll]*R[cells[ii][ll]][1];
01264                                       a3[kk]=a3[kk]+Q[kk][ll]*R[cells[ii][ll]][2];
01265                                    }
01266                                 }
01267                                 det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
01268                                 if(gqmin>det) gqmin=det;
01269                                     if((i==nn)&&(j==l)&&(k==m)) sigma=(double)1/27;
01270                                 if(((i==nn)&&(j==l)&&(k!=m))||((i==nn)&&(j!=l)&&(k==m))||((i!=nn)&&(j==l)&&(k==m))) sigma=(double)1/(27*2);
01271                                 if(((i==nn)&&(j!=l)&&(k!=m))||((i!=nn)&&(j!=l)&&(k==m))||((i!=nn)&&(j==l)&&(k!=m))) sigma=(double)1/(27*4);
01272                                 if((i!=nn)&&(j!=l)&&(k!=m)) sigma=(double)1/(27*8);
01273                                 v+=sigma*det;
01274                                 vcell+=sigma*det;
01275                             }
01276                          }
01277                       }
01278                   }
01279               }
01280           }
01281           if(vmin>vcell) vmin=vcell;
01282        } else {//quad tetr
01283           vcell=0;
01284           for(i=0;i<4;i++){
01285               for(j=0;j<4;j++){
01286                   for(k=0;k<4;k++){
01287                     switch (i)
01288                     { case 0 : K[0]=0; K[1]=0; K[2]=0; break;
01289                       case 1 : K[0]=1; K[1]=0; K[2]=0; break;
01290                       case 2 : K[0]=1.0/2; K[1]=1; K[2]=0; break;
01291                       case 3 : K[0]=1.0/2; K[1]=1.0/3; K[2]=1; break; }
01292                     switch (j)
01293                     { case 0 : K[3]=0; K[4]=0; K[5]=0; break;
01294                       case 1 : K[3]=1; K[4]=0; K[5]=0; break;
01295                       case 2 : K[3]=1.0/2; K[4]=1; K[5]=0; break;
01296                       case 3 : K[3]=1.0/2; K[4]=1.0/3; K[5]=1; break; }
01297                     switch (k)
01298                     { case 0 : K[6]=0; K[7]=0; K[8]=0; break;
01299                       case 1 : K[6]=1; K[7]=0; K[8]=0; break;
01300                       case 2 : K[6]=1.0/2; K[7]=1; K[8]=0; break;
01301                       case 3 : K[6]=1.0/2; K[7]=1.0/3; K[8]=1; break; }
01302                     basisA(3,Q,10,K,H[ii],me);
01303                     for(kk=0;kk<3;kk++){a1[kk]=0; a2[kk]=0; a3[kk]=0;
01304                        for(ll=0;ll<10;ll++){
01305                           a1[kk]=a1[kk]+Q[kk][ll]*R[cells[ii][ll]][0];
01306                           a2[kk]=a2[kk]+Q[kk][ll]*R[cells[ii][ll]][1];
01307                           a3[kk]=a3[kk]+Q[kk][ll]*R[cells[ii][ll]][2];
01308                        }
01309                     }
01310                     det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
01311                     if(gqmin>det) gqmin=det;
01312                     if((i==j)&&(j==k)) sigma=1.0/120; else
01313                                                 if((i==j)||(j==k)||(i==k)) sigma=1.0/360; else sigma=1.0/720;
01314                     v+=sigma*det;   vcell+=sigma*det;
01315                           }
01316                           }
01317                   }
01318                   if(vmin>vcell) vmin=vcell;
01319        }
01320     }
01321 }
01322 
01323 
01324 (*vol)=v/(double)ncells;
01325 (*Vmin)=vmin;
01326 
01327   for(i=0;i<n;i++) free(Q[i]);
01328   free(Q);
01329 
01330 return gqmin;
01331 }

int libMesh::VariationalMeshSmoother::pcg_ic0	(	int	n,
		LPINT	ia,
		LPINT	ja,
		LPDOUBLE	a,
		LPDOUBLE	u,
		LPDOUBLE	x,
		LPDOUBLE	b,
		LPDOUBLE	r,
		LPDOUBLE	p,
		LPDOUBLE	z,
		double	eps,
		int	maxite,
		int	msglev,
		FILE *	sout
	)			`[private]`

Solve the SPD linear system by PCG method

Definition at line 2526 of file mesh_smoother_vsmoother.C.

Referenced by solver().

02528 {
02529     int i, ii, j, k;
02530     double beta, pap, rhr = 0.0, rhr0 = 0.0, rhrold = 0.0, alpha, rr0 = 0.0, rr;
02531 
02532     for(i=0;i<=maxite;i++){
02533            if(i==0) for(k=0;k<n;k++) p[k]=x[k];
02534            //z=Ap
02535            for(ii=0;ii<n;ii++) z[ii]=0.0;
02536        for(ii=0;ii<n;ii++){
02537                    z[ii]+=a[ia[ii]]*p[ii];
02538                    for(j=ia[ii]+1;j<ia[ii+1];j++){
02539                        z[ii]+=a[j]*p[ja[j]-1];
02540                        z[ja[j]-1]+=a[j]*p[ii];
02541                            }
02542            }
02543            if(i==0) for(k=0;k<n;k++) r[k]=b[k]-z[k]; //r(0) = b - Ax(0)
02544            if(i>0){
02545                pap=0.0;
02546                for(k=0;k<n;k++) pap+=p[k]*z[k];
02547                alpha=rhr/pap;
02548            for(k=0;k<n;k++){
02549                        x[k]+=p[k]*alpha;
02550                        r[k]-=z[k]*alpha;
02551                    }
02552                rhrold = rhr;
02553            }
02554        // z = H * r
02555        for(ii=0;ii<n;ii++) z[ii]=r[ii]*u[ii];
02556            if(i==0) for(k=0;k<n;k++) p[k]=z[k];// p(0) = Hr(0)
02557            rhr=0.0;
02558            rr=0.0;
02559            for(k=0;k<n;k++){
02560                rhr+=r[k]*z[k];
02561                rr+=r[k]*r[k];
02562            }
02563            if(i==0){
02564             rhr0 = rhr;
02565             rr0 = rr;
02566            }
02567            if(msglev>0) fprintf(sout, "%d )  rHr =%g  rr =%g \n", i, sqrt(rhr/rhr0), sqrt(rr/rr0));
02568 
02569            if(rr<=eps*eps*rr0) return i;
02570 
02571            if(i>=maxite) return i;
02572 
02573            if(i>0){
02574               beta=rhr/rhrold;
02575               for(k=0;k<n;k++) p[k]=z[k]+p[k]*beta;
02576            }
02577         }
02578 
02579     return i;
02580 }

int libMesh::VariationalMeshSmoother::pcg_par_check	(	int	n,
		LPINT	ia,
		LPINT	ja,
		LPDOUBLE	a,
		double	eps,
		int	maxite,
		int	msglev,
		FILE *	sout
	)			`[private]`

Input parameter check

Definition at line 2462 of file mesh_smoother_vsmoother.C.

Referenced by solver().

02463 {
02464     int i, j, jj, k;
02465 
02466     if (n<=0){
02467             if (msglev>0) fprintf(sout, "sol_pcg: incorrect input parameter: n =%d \n",n);
02468             return -1;
02469     }
02470     if (ia[0]!=0) {
02471             if (msglev > 0) fprintf(sout, "sol_pcg: incorrect input parameter: ia(1) =%d \n", ia[0]);
02472             return -2;
02473     }
02474     for (i=0;i<n;i++) {
02475             if (ia[i+1]<ia[i]) {
02476             if (msglev>=1)fprintf(sout, "sol_pcg: incorrect input parameter: i =%d ; ia(i) =%d  must be <= a(i+1) =%d \n",
02477                         (i+1), ia[i], ia[i+1]);
02478             return -2;
02479                 }
02480         }
02481     for (i=0;i<n;i++) {
02482             if (ja[ia[i]]!=(i+1)) {
02483             if (msglev>=1) fprintf(sout, "sol_pcg: incorrect input parameter: i =%d ; ia(i) =%d ; absence of diagonal entry; k =%d ; the value ja(k) =%d  must be equal to i\n",
02484                         (i+1), ia[i], ia[i]+1, ja[ia[i]]);
02485                 return -3;
02486                 }
02487             jj = 0;
02488             for (k=ia[i];k<ia[i+1];k++) {
02489                 j=ja[k];
02490                 if ((j<(i+1))||(j>n)) {
02491                     if (msglev>=1) fprintf(sout, "sol_pcg: incorrect input parameter: i =%d ; ia(i) =%d ; a(i+1) =%d ; k =%d ; the value ja(k) =%d  is out of range\n",
02492                             (i+1), ia[i], ia[i+1], (k+1), ja[k]);
02493                     return -3;
02494                         }
02495                 if (j<=jj) {
02496                        if (msglev >= 1) fprintf(sout, "sol_pcg: incorrect input parameter: i =%d ; ia(i) =%d ; a(i+1) =%d ; k =%d ; the value ja(k) =%d  must be less than ja(k+1) =%d \n",
02497                                                (i+1), ia[i], ia[i+1], (k+1), ja[k], ja[k+1]);
02498                     return -3;
02499                         }
02500                         jj=j;
02501                 }
02502         }
02503 
02504     for (i=0;i<n;i++) {
02505           if (a[ia[i]]<=0.0e0) {
02506             if (msglev>=1) fprintf(sout, "sol_pcg: incorrect input parameter: i =%d ; ia(i) =%d ; the diagonal entry a(ia(i)) =%g  must be > 0\n",
02507                         (i+1), ia[i], a[ia[i]]);
02508             return -4;
02509           }
02510         }
02511 
02512     if (eps<0.0e0) {
02513             if (msglev>0) fprintf(sout, "sol_pcg: incorrect input parameter: eps =%g \n",eps);
02514             return -7;
02515     }
02516     if (maxite<1) {
02517             if (msglev>0) fprintf(sout, "sol_pcg: incorrect input parameter: maxite =%d \n", maxite);
02518             return -8;
02519     }
02520     return 0;
02521 }

int libMesh::VariationalMeshSmoother::read_adp	(	LPDOUBLE	afun,
		char *	adap,
		FILE *	sout
	)			`[private]`

int libMesh::VariationalMeshSmoother::readgr	(	int	n,
		int	N,
		LPLPDOUBLE	R,
		LPINT	mask,
		int	ncells,
		LPLPINT	cells,
		LPINT	mcells,
		int	nedges,
		LPINT	edges,
		LPINT	hnodes,
		FILE *	sout
	)			`[private]`

reading grid from input file

Definition at line 279 of file mesh_smoother_vsmoother.C.

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

Referenced by metr_data_gen(), and smooth().

00281 {
00282   libMesh::out<<"Sarting readgr"<<std::endl;
00283   int i;
00284  //add error messages where format can be inconsistent
00285 
00286   //Find the boundary nodes
00287   std::vector<bool> on_boundary;
00288   MeshTools::find_boundary_nodes(_mesh, on_boundary);
00289 
00290   //Grab node coordinates and set mask
00291   {
00292     MeshBase::const_node_iterator       it  = _mesh.nodes_begin();
00293     const MeshBase::const_node_iterator end = _mesh.nodes_end();
00294 
00295     //Only compute the node to elem map once
00296     std::vector<std::vector<const Elem*> > nodes_to_elem_map;
00297     MeshTools::build_nodes_to_elem_map(_mesh, nodes_to_elem_map);
00298 
00299     i=0;
00300     for (; it != end; ++it)
00301     {
00302       //For each node grab its X Y [Z] coordinates
00303       for(unsigned int j=0;j<_dim;j++)
00304         R[i][j]=(*(*it))(j);
00305 
00306       //Set the Proper Mask
00307       //Internal nodes are 0
00308       //Immovable boundary nodes are 1
00309       //Movable boundary nodes are 2
00310       if(on_boundary[i])
00311       {
00312         //Only look for sliding edge nodes in 2D
00313         if(_dim == 2)
00314         {
00315           //Find all the nodal neighbors... that is the nodes directly connected
00316           //to this node through one edge
00317           std::vector<const Node*> neighbors;
00318           MeshTools::find_nodal_neighbors(_mesh, *(*it), nodes_to_elem_map, neighbors);
00319 
00320           std::vector<const Node*>::const_iterator ne = neighbors.begin();
00321           std::vector<const Node*>::const_iterator ne_end = neighbors.end();
00322 
00323           //Grab the x,y coordinates
00324           Real x = (*(*it))(0);
00325           Real y = (*(*it))(1);
00326 
00327           //Theta will represent the atan2 angle (meaning with the proper quadrant in mind)
00328           //of the neighbor node in a system where the current node is at the origin
00329           Real theta = 0;
00330           std::vector<Real> thetas;
00331 
00332           //Calculate the thetas
00333           for(; ne != ne_end; ne++)
00334           {
00335             //Note that the x and y values of this node are subtracted off
00336             //this centers the system around this node
00337             theta = atan2((*(*ne))(1)-y,(*(*ne))(0)-x);
00338 
00339             //Save it for later
00340             thetas.push_back(theta);
00341           }
00342 
00343           //Assume the node is immovable... then prove otherwise
00344           mask[i]=1;
00345 
00346           //Search through neighbor nodes looking for two that form a straight line with this node
00347           for(unsigned int a=0;a<thetas.size()-1;a++)
00348           {
00349             //Only try each pairing once
00350             for(unsigned int b=a+1;b<thetas.size();b++)
00351             {
00352               //Find if the two neighbor nodes angles are 180 degrees (pi) off of eachother (withing a tolerance)
00353               //In order to make this a true movable boundary node... the two that forma  straight line with
00354               //it must also be on the boundary
00355               if(on_boundary[neighbors[a]->id()] &&
00356                 on_boundary[neighbors[b]->id()] &&
00357                 (
00358                   (fabs(thetas[a]-(thetas[b] + (libMesh::pi))) < .001) ||
00359                   (fabs(thetas[a]-(thetas[b] - (libMesh::pi))) < .001)
00360                 )
00361                 )
00362               {
00363                 //if((*(*it))(1) > 0.25 || (*(*it))(0) > .7 || (*(*it))(0) < .01)
00364                   mask[i]=2;
00365               }
00366 
00367             }
00368           }
00369         }
00370         else //In 3D set all boundary nodes to be fixed
00371           mask[i]=1;
00372       }
00373       else
00374         mask[i]=0;  //Internal Node
00375 
00376       //libMesh::out<<"Node: "<<i<<"  Mask: "<<mask[i]<<std::endl;
00377       i++;
00378     }
00379   }
00380 
00381   // Grab the connectivity
00382   // FIXME: Generalize this!
00383   {
00384     MeshBase::const_element_iterator it  = _mesh.active_elements_begin();
00385     const MeshBase::const_element_iterator end = _mesh.active_elements_end();
00386 
00387     i=0;
00388     for (; it != end; ++it)
00389     {
00390       //Keep track of the number of nodes
00391       //there must be 6 for 2D
00392       //10 for 3D
00393       //If there are actually less than that -1 must be used
00394       //to fill out the rest
00395       int num=0;
00396 /*
00397       int num_necessary = 0;
00398 
00399       if (_dim == 2)
00400         num_necessary = 6;
00401       else
00402         num_necessary = 10;
00403 */
00404 
00405       if(_dim == 2)
00406       {
00407         switch((*it)->n_vertices())
00408         {
00409           //Grab nodes that do exist
00410           case 3:  //Tri
00411             for(unsigned int k=0;k<(*it)->n_vertices();k++)
00412               cells[i][k]=(*it)->node(k);
00413             num=(*it)->n_vertices();
00414             break;
00415           case 4:  //Quad 4
00416             cells[i][0]=(*it)->node(0);
00417             cells[i][1]=(*it)->node(1);
00418             cells[i][2]=(*it)->node(3); //Note that 2 and 3 are switched!
00419             cells[i][3]=(*it)->node(2);
00420             num=4;
00421             break;
00422           default:
00423             libmesh_assert(false);
00424         }
00425       }
00426       else
00427       {
00428         //Grab nodes that do exist
00429         switch((*it)->n_vertices())
00430         {
00431           case 4:  //Tet 4
00432             for(unsigned int k=0;k<(*it)->n_vertices();k++)
00433               cells[i][k]=(*it)->node(k);
00434             num=(*it)->n_vertices();
00435             break;
00436           case 8:  //Hex 8f
00437             cells[i][0]=(*it)->node(0);
00438             cells[i][1]=(*it)->node(1);
00439             cells[i][2]=(*it)->node(3); //Note that 2 and 3 are switched!
00440             cells[i][3]=(*it)->node(2);
00441 
00442             cells[i][4]=(*it)->node(4);
00443             cells[i][5]=(*it)->node(5);
00444             cells[i][6]=(*it)->node(7); //Note that 6 and 7 are switched!
00445             cells[i][7]=(*it)->node(6);
00446             num=8;
00447           default:
00448             libmesh_assert(false);
00449         }
00450       }
00451 
00452       //Fill in the rest with -1
00453       for(int j=num;j<3*n+n%2;j++)
00454         cells[i][j]=-1;
00455 
00456       //Mask it with 0 to state that this is an active element
00457       //FIXME: Could be something other than zero
00458       mcells[i]=0;
00459 
00460       i++;
00461     }
00462   }
00463 
00464   //Grab hanging node connectivity
00465   {
00466     std::map<dof_id_type, std::vector<dof_id_type> >::iterator it = _hanging_nodes.begin();
00467     std::map<dof_id_type, std::vector<dof_id_type> >::iterator end = _hanging_nodes.end();
00468 
00469     for(i=0;it!=end;it++)
00470     {
00471 
00472       libMesh::out<<"Parent 1: "<<(it->second)[0]<<std::endl;
00473       libMesh::out<<"Parent 2: "<<(it->second)[1]<<std::endl;
00474       libMesh::out<<"Hanging Node: "<<it->first<<std::endl<<std::endl;
00475 
00476 
00477       //First Parent
00478       edges[2*i]=(it->second)[1];
00479 
00480       //Second Parent
00481       edges[2*i+1]=(it->second)[0];
00482 
00483       //Hanging Node
00484       hnodes[i]=it->first;
00485 
00486       i++;
00487     }
00488   }
00489   libMesh::out<<"Finished readgr"<<std::endl;
00490 
00491 return 0;
00492 }

int libMesh::VariationalMeshSmoother::readmetr	(	char *	name,
		LPLPLPDOUBLE	H,
		int	ncells,
		int	n,
		FILE *	sout
	)			`[private]`

Read Metrics

Definition at line 498 of file mesh_smoother_vsmoother.C.

Referenced by smooth().

00499 {
00500 int i, j, k, scanned;
00501 double d;
00502 FILE *stream;
00503 
00504    stream=fopen(name,"r");
00505    for(i=0;i<ncells;i++)
00506           for(j=0;j<n;j++){
00507           for(k=0;k<n;k++){scanned = fscanf(stream,"%le ",&d);
00508             libmesh_assert_not_equal_to (scanned, EOF);
00509             H[i][j][k]=d;}
00510                   scanned = fscanf(stream,"\n");
00511                   libmesh_assert_not_equal_to (scanned, EOF);
00512           }
00513 
00514    fclose(stream);
00515 
00516 return 0;
00517 }

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

References _adaptive_func, _dim, _dist_norm, _generate_data, _hanging_nodes, _maxiter, libMesh::MeshSmoother::_mesh, _metric, _miniter, _miniterBC, _theta, alloc_d_n1(), alloc_d_n1_n2(), alloc_d_n1_n2_n3(), alloc_i_n1(), alloc_i_n1_n2(), 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().

00046 {
00047   int n, me, gr, adp, miniter, miniterBC, maxiter, N, ncells, nedges, s, vms_err=0;
00048   double theta;
00049   char grid[40], metr[40], grid_old[40], adap[40];
00050 //  FILE *stream;
00051   LPLPLPDOUBLE H;
00052   LPLPDOUBLE R;
00053   LPLPINT cells;
00054   LPINT mask, edges, mcells, hnodes;
00055   int iter[4],i,j;
00056   clock_t ticks1, ticks2;
00057 
00058   FILE *sout;
00059   sout=fopen("smoother.out","wr");//stdout;
00060 
00061   n=_dim;
00062   me=_metric;
00063   gr=_generate_data?0:1;
00064   adp=_adaptive_func;
00065   theta=_theta;
00066   miniter=_miniter;
00067   maxiter=_maxiter;
00068   miniterBC=_miniterBC;
00069 
00070   if((gr==0)&&(me>1))
00071   {
00072     for(i=0;i<40;i++) grid_old[i]=grid[i];
00073       metr_data_gen(grid, metr, n, me, sout); //generate metric from initial mesh (me=2,3)
00074   }
00075 
00076   s=_dim;
00077   N=_mesh.n_nodes();
00078   ncells=_mesh.n_active_elem();
00079 
00080   //I wish I could do this in readgr... but the way she allocs
00081   //memory we need to do this here... or pass a lot of things around
00082   MeshTools::find_hanging_nodes_and_parents(_mesh,_hanging_nodes);
00083 
00084   nedges=_hanging_nodes.size();
00085   //nedges=0;
00086   if(s!=n)
00087     {
00088       fprintf(sout,"Error: dim in input file is inconsistent with dim in .cfg \n");
00089       fclose(sout);
00090       return _dist_norm;
00091     }
00092 
00093   mask=alloc_i_n1(N);
00094   edges=alloc_i_n1(2*nedges);
00095   mcells=alloc_i_n1(ncells);
00096   hnodes=alloc_i_n1(nedges);
00097   H=alloc_d_n1_n2_n3(ncells,n,n);
00098   R=alloc_d_n1_n2(N,n);
00099   cells=alloc_i_n1_n2(ncells,3*n+n%2);
00100   for(i=0;i<ncells;i++){
00101      cells[i]=alloc_i_n1(3*n+n%2);
00102          if(me>1){
00103                  H[i]=alloc_d_n1_n2(n,n);
00104                  for(j=0;j<n;j++) H[i][j]=alloc_d_n1(n);
00105          }
00106   }
00107   for(i=0;i<N;i++) R[i]=alloc_d_n1(n);
00108 
00109   /*----------initial grid---------*/
00110   vms_err=readgr(n,N,R,mask,ncells,cells,mcells,nedges,edges,hnodes,sout);
00111   if(vms_err<0)
00112     {
00113       fprintf(sout,"Error reading input mesh file\n");
00114       fclose(sout);
00115       return _dist_norm;
00116     }
00117   if(me>1) vms_err=readmetr(metr,H,ncells,n,sout);
00118   if(vms_err<0)
00119     { 
00120       fprintf(sout,"Error reading metric file\n");
00121       fclose(sout);
00122       return _dist_norm;
00123     }
00124 
00125   iter[0]=miniter; iter[1]=maxiter; iter[2]=miniterBC;
00126 
00127 
00128   /*---------grid optimization--------*/
00129   fprintf(sout,"Starting Grid Optimization \n");
00130   ticks1=clock();
00131   full_smooth(n,N,R,mask,ncells,cells,mcells,nedges,edges,hnodes,theta,iter,me,H,adp,adap,gr,sout);
00132   ticks2=clock();
00133   fprintf(sout,"full_smooth took (%d-%d)/%ld = %ld seconds \n",
00134           (int)ticks2,(int)ticks1,(long int)CLOCKS_PER_SEC,(long int)(ticks2-ticks1)/CLOCKS_PER_SEC);
00135 
00136   /*---------save result---------*/
00137   fprintf(sout,"Saving Result \n");
00138   writegr(n,N,R,mask,ncells,cells,mcells,nedges,edges,hnodes,"fin_grid.dat", 4-me+3*gr, grid_old, sout);
00139 
00140   for(i=0;i<N;i++) free(R[i]);
00141   for(i=0;i<ncells;i++){
00142           if(me>1){
00143                   for(j=0;j<n;j++) free(H[i][j]);
00144                   free(H[i]);
00145           }
00146      free(cells[i]);
00147   }
00148 
00149   free(cells);
00150   free(H);
00151   free(R);
00152   free(mask);
00153   free(edges);
00154   free(mcells);
00155   free(hnodes);
00156   fclose(sout);
00157   libmesh_assert_greater (_dist_norm, 0);
00158 
00159   return _dist_norm;
00160 }

virtual void libMesh::VariationalMeshSmoother::smooth ( ) [inline, 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 166 of file mesh_smoother_vsmoother.h.

References _distance, and smooth().

Referenced by smooth().

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

int libMesh::VariationalMeshSmoother::solver	(	int	n,
		LPINT	ia,
		LPINT	ja,
		LPDOUBLE	a,
		LPDOUBLE	x,
		LPDOUBLE	b,
		double	eps,
		int	maxite,
		int	msglev,
		FILE *	sout
	)			`[private]`

Solve Symmetric Positive Definite (SPD) linear system by Conjugate Gradient (CG) method preconditioned by Point Jacobi (diagonal scaling)

Definition at line 2418 of file mesh_smoother_vsmoother.C.

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

Referenced by minJ().

02419 {
02420     int info, i;
02421     LPDOUBLE u, r, p, z;
02422 
02423     info = 0;
02424 
02425     fprintf(sout, "Beginning Solve %d\n", n);
02426 
02427 
02428     // Check parameters
02429     info=pcg_par_check(n, ia, ja, a, eps, maxite, msglev, sout);
02430     if (info != 0)
02431         return info;
02432 
02433     // Allocate working arrays
02434     u=alloc_d_n1(n);
02435     r=alloc_d_n1(n);
02436     p=alloc_d_n1(n);
02437     z=alloc_d_n1(n);
02438 
02439     // PJ preconditioner construction
02440         for (i=0;i<n;i++) u[i]=1.0/a[ia[i]];
02441 
02442     // PCG iterations
02443     info=pcg_ic0(n, ia, ja, a, u, x, b, r, p, z, eps, maxite, msglev, sout);
02444 
02445     // Free working arrays
02446     free(u);
02447         free(r);
02448         free(p);
02449         free(z);
02450 
02451     //Mat sparse_global;
02452     //MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,n,n,ia,ja,a,&sparse_global);
02453 
02454     fprintf(sout, "Finished Solve %d\n",n);
02455 
02456     return info;
02457 }

double libMesh::VariationalMeshSmoother::vertex	(	int	n,
		LPLPLPDOUBLE	W,
		LPLPDOUBLE	F,
		LPLPDOUBLE	R,
		LPINT	cell_in,
		double	epsilon,
		double	w,
		int	nvert,
		LPDOUBLE	K,
		LPLPDOUBLE	H,
		int	me,
		double	vol,
		int	f,
		double *	qmin,
		int	adp,
		LPDOUBLE	g,
		double	sigma,
		FILE *	sout
	)			`[private]`

Computes local matrics W and local rhs F on one basis

Definition at line 2169 of file mesh_smoother_vsmoother.C.

References alloc_d_n1(), alloc_d_n1_n2(), alloc_d_n1_n2_n3(), basisA(), jac2(), jac3(), and std::pow().

Referenced by localP().

02173 {
02174   LPLPLPDOUBLE P = NULL, d2phi = NULL;
02175   LPLPDOUBLE  gpr = NULL, dphi = NULL, dfe = NULL;
02176   LPLPDOUBLE Q;
02177   double a1[3], a2[3], a3[3], av1[3], av2[3], av3[3];
02178   int i,j,k,i1;
02179   double tr=0.0, det=0.0, dchi, chi, fet, phit=0.0, G, con=100.0;
02180   Q = alloc_d_n1_n2(3, nvert);
02181   for(i=0;i<n;i++) Q[i]=alloc_d_n1(nvert);
02182 if(f==0){
02183   gpr = alloc_d_n1_n2(3, nvert);
02184   dphi=alloc_d_n1_n2(3,3);
02185   dfe=alloc_d_n1_n2(3,3);
02186   for(i=0;i<n;i++){ gpr[i]=alloc_d_n1(nvert);
02187       dphi[i]=alloc_d_n1(3);
02188       dfe[i]=alloc_d_n1(3);}
02189   P = alloc_d_n1_n2_n3(3,3,3);
02190   d2phi = alloc_d_n1_n2_n3(3,3,3);
02191   for(i=0;i<n;i++){ P[i]=alloc_d_n1_n2(3,3);
02192      d2phi[i]=alloc_d_n1_n2(3,3);
02193      for(j=0;j<n;j++){ P[i][j]=alloc_d_n1(3);
02194         d2phi[i][j]=alloc_d_n1(3);}
02195   }
02196 }
02197 
02198   /*--------hessian, function, gradient-----------------------*/
02199   basisA(n,Q,nvert,K,H,me);
02200   for(i=0;i<n;i++)
02201   {
02202     a1[i]=0;
02203     a2[i]=0;
02204     a3[i]=0;
02205     for(j=0;j<nvert;j++)
02206     {
02207       a1[i]+=Q[i][j]*R[cell_in[j]][0];
02208       a2[i]+=Q[i][j]*R[cell_in[j]][1];
02209       if(n==3)
02210         a3[i]+=Q[i][j]*R[cell_in[j]][2];
02211     }
02212   }
02213   //account for adaptation
02214   if(adp<2)
02215     G=g[0];
02216   else
02217   {
02218     G=1.0;
02219     for(i=0;i<n;i++)
02220     {
02221       a1[i]*=sqrt(g[0]);
02222       a2[i]*=sqrt(g[1]);
02223       if(n==3)
02224         a3[i]*=sqrt(g[2]);
02225     }
02226   }
02227 
02228   if(n==2)
02229   {
02230     av1[0]= a2[1];  av1[1]=-a2[0];
02231     av2[0]=-a1[1];  av2[1]= a1[0];
02232     det=jac2(a1[0],a1[1],a2[0],a2[1]);
02233     tr=0;
02234     for(i=0;i<2;i++) tr=tr+a1[i]*a1[i]/2+a2[i]*a2[i]/2;
02235       phit=(1-w)*tr+w*(vol/G+det*det*G/vol)/(double)2;
02236   }
02237 
02238   if(n==3)
02239   {
02240       av1[0]= a2[1]*a3[2]-a2[2]*a3[1];  av1[1]=a2[2]*a3[0]-a2[0]*a3[2]; av1[2]=a2[0]*a3[1]-a2[1]*a3[0];
02241       av2[0]= a3[1]*a1[2]-a3[2]*a1[1];  av2[1]=a3[2]*a1[0]-a3[0]*a1[2]; av2[2]=a3[0]*a1[1]-a3[1]*a1[0];
02242       av3[0]= a1[1]*a2[2]-a1[2]*a2[1];  av3[1]=a1[2]*a2[0]-a1[0]*a2[2]; av3[2]=a1[0]*a2[1]-a1[1]*a2[0];
02243           det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
02244       tr=0;
02245       for(i=0;i<3;i++) tr=tr+a1[i]*a1[i]/3+a2[i]*a2[i]/3+a3[i]*a3[i]/3;
02246       phit=(1-w)*pow(tr,1.5)+w*(vol/G+det*det*G/vol)/(double)2;
02247   }
02248 
02249   dchi=0.5+det*0.5/sqrt(epsilon*epsilon+det*det);
02250   chi=0.5*(det+sqrt(epsilon*epsilon+det*det));
02251   if(chi!=0) fet=phit/chi; else fet=1e32; //function is computed
02252   (*qmin)=det*G;
02253 
02254   //gradient and Hessian
02255   if(f==0)
02256   {
02257     if(n==2)
02258     {
02259       for(i=0;i<2;i++)
02260       {
02261         dphi[0][i]=(1-w)*a1[i]+w*G*det*av1[i]/vol;
02262         dphi[1][i]=(1-w)*a2[i]+w*G*det*av2[i]/vol;
02263       }
02264 
02265       for(i=0;i<2;i++)
02266       {
02267         for(j=0;j<2;j++)
02268         {
02269           d2phi[0][i][j]=w*G*av1[i]*av1[j]/vol;
02270           d2phi[1][i][j]=w*G*av2[i]*av2[j]/vol;
02271 
02272           if(i==j) for(k=0;k<2;k++)
02273           {
02274             d2phi[k][i][j]+=(1-w);
02275           }
02276         }
02277       }
02278 
02279       for(i=0;i<2;i++)
02280       {
02281         dfe[0][i]=(dphi[0][i]-fet*dchi*av1[i])/chi;
02282         dfe[1][i]=(dphi[1][i]-fet*dchi*av2[i])/chi;
02283       }
02284 
02285       for(i=0;i<2;i++)
02286       {
02287         for(j=0;j<2;j++)
02288         {
02289           P[0][i][j]=(d2phi[0][i][j]-dchi*(av1[i]*dfe[0][j]+dfe[0][i]*av1[j]))/chi;
02290           P[1][i][j]=(d2phi[1][i][j]-dchi*(av2[i]*dfe[1][j]+dfe[1][i]*av2[j]))/chi;
02291         }
02292       }
02293     }
02294 
02295     if(n==3)
02296     {
02297       for(i=0;i<3;i++)
02298       {
02299         dphi[0][i]=(1-w)*sqrt(tr)*a1[i]+w*G*det*av1[i]/vol;
02300         dphi[1][i]=(1-w)*sqrt(tr)*a2[i]+w*G*det*av2[i]/vol;
02301         dphi[2][i]=(1-w)*sqrt(tr)*a3[i]+w*G*det*av3[i]/vol;
02302       }
02303 
02304       for(i=0;i<3;i++)
02305       {
02306         if(tr!=0)
02307         {
02308           d2phi[0][i][i]=(1-w)/((double)3*sqrt(tr))*a1[i]*a1[i]+w*G*av1[i]*av1[i]/vol;
02309           d2phi[1][i][i]=(1-w)/((double)3*sqrt(tr))*a2[i]*a2[i]+w*G*av2[i]*av2[i]/vol;
02310           d2phi[2][i][i]=(1-w)/((double)3*sqrt(tr))*a3[i]*a3[i]+w*G*av3[i]*av3[i]/vol;
02311         }
02312         else
02313         {
02314           for(k=0;k<3;k++)
02315             d2phi[k][i][i]=0;
02316         }
02317 
02318         for(k=0;k<3;k++)
02319           d2phi[k][i][i]+=(1-w)*sqrt(tr);
02320       }
02321 
02322       for(i=0;i<3;i++)
02323       {
02324         dfe[0][i]=(dphi[0][i]-fet*dchi*av1[i]+con*epsilon*a1[i])/chi;
02325         dfe[1][i]=(dphi[1][i]-fet*dchi*av2[i]+con*epsilon*a2[i])/chi;
02326         dfe[2][i]=(dphi[2][i]-fet*dchi*av3[i]+con*epsilon*a3[i])/chi;
02327       }
02328 
02329       for(i=0;i<3;i++)
02330       {
02331         P[0][i][i]=(d2phi[0][i][i]-dchi*(av1[i]*dfe[0][i]+dfe[0][i]*av1[i])+con*epsilon)/chi;
02332         P[1][i][i]=(d2phi[1][i][i]-dchi*(av2[i]*dfe[1][i]+dfe[1][i]*av2[i])+con*epsilon)/chi;
02333         P[2][i][i]=(d2phi[2][i][i]-dchi*(av3[i]*dfe[2][i]+dfe[2][i]*av3[i])+con*epsilon)/chi;
02334       }
02335 
02336       for(k=0;k<3;k++)
02337       {
02338         for(i=0;i<3;i++)
02339         {
02340           for(j=0;j<3;j++)
02341           {
02342             if(i!=j)
02343               P[k][i][j]=0;
02344           }
02345         }
02346       }
02347     }
02348 
02349     /*--------matrix W, right side F---------------------*/
02350     for(i1=0;i1<n;i1++)
02351     {
02352       for(i=0;i<n;i++)
02353       {
02354         for(j=0;j<nvert;j++)
02355         {
02356           gpr[i][j]=0;
02357           for(k=0;k<n;k++) gpr[i][j]=gpr[i][j]+P[i1][i][k]*Q[k][j];
02358         }
02359       }
02360 
02361       for(i=0;i<nvert;i++)
02362       {
02363         for(j=0;j<nvert;j++)
02364         {
02365           for(k=0;k<n;k++)
02366           {
02367             W[i1][i][j]=W[i1][i][j]+Q[k][i]*gpr[k][j]*sigma;
02368           }
02369         }
02370       }
02371 
02372       for(i=0;i<nvert;i++)
02373       {
02374         for(k=0;k<2;k++)
02375           F[i1][i]=F[i1][i]+Q[k][i]*dfe[i1][k]*sigma;
02376       }
02377     }
02378   }
02379 
02380   /*----------------------------------------*/
02381   for(i=0;i<n;i++)
02382     free(Q[i]);
02383   free(Q);
02384   if(f==0)
02385   {
02386     for(i=0;i<n;i++)
02387     {
02388       free(gpr[i]);
02389       free(dphi[i]);
02390       free(dfe[i]);
02391     }
02392     free(gpr);
02393     free(dphi);
02394     free(dfe);
02395 
02396     for(i=0;i<n;i++)
02397     {
02398       for(j=0;j<n;j++)
02399       {
02400         free(P[i][j]);
02401         free(d2phi[i][j]);
02402       }
02403       free(P[i]);
02404       free(d2phi[i]);
02405     }
02406     free(P);
02407     free(d2phi);
02408   }
02409 
02410   return (fet*sigma);
02411 }

int libMesh::VariationalMeshSmoother::writegr	(	int	n,
		int	N,
		LPLPDOUBLE	R,
		LPINT	mask,
		int	ncells,
		LPLPINT	cells,
		LPINT	mcells,
		int	nedges,
		LPINT	edges,
		LPINT	hnodes,
		const char	grid[],
		int	me,
		const char	grid_old[],
		FILE *	sout
	)			`[private]`

save grid

Definition at line 167 of file mesh_smoother_vsmoother.C.

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

Referenced by smooth().

00169 {
00170   libMesh::out<<"Starting writegr"<<std::endl;
00171   int i;
00172 
00173   //Adjust nodal coordinates to new positions
00174   {
00175     MeshBase::node_iterator       it  = _mesh.nodes_begin();
00176     const MeshBase::node_iterator end = _mesh.nodes_end();
00177 
00178     libmesh_assert_equal_to (_dist_norm, 0.0);
00179     _dist_norm=0;
00180     i=0;
00181     for (; it != end; ++it)
00182     {
00183       double total_dist=0.0;
00184 
00185       //For each node set its X Y [Z] coordinates
00186       for(unsigned int j=0;j<_dim;j++)
00187       {
00188         double distance = R[i][j]-(*(*it))(j);
00189 
00190         //Save the squares of the distance
00191         total_dist+=Utility::pow<2>(distance);
00192 
00193         (*(*it))(j)=(*(*it))(j)+(distance*_percent_to_move);
00194       }
00195 
00196       libmesh_assert_greater_equal (total_dist, 0.0);
00197 
00198       //Add the distance this node moved to the global distance
00199       _dist_norm+=total_dist;
00200 
00201       i++;
00202     }
00203 
00204     //Relative "error"
00205     _dist_norm=sqrt(_dist_norm/_mesh.n_nodes());
00206   }
00207   /*
00208 int scanned;
00209 if(me>=3){
00210         fprintf(stream,"%d \n%d \n%d \n%d \n",n,N,ncells,nedges);
00211 
00212     for(i=0;i<N;i++){//node coordinates
00213             for(unsigned int j=0;j<n;j++) fprintf(stream,"%e ",R[i][j]);
00214         fprintf(stream,"%d \n",mask[i]);
00215     }
00216     for(i=0;i<ncells;i++){//cell connectivity
00217             int nvert=0;
00218             while(cells[i][nvert]>=0){
00219                fprintf(stream,"%d ",cells[i][nvert]);
00220                nvert++;
00221             }
00222             for(unsigned int j=nvert;j<3*n+n%2;j++) fprintf(stream,"-1 ");
00223             fprintf(stream,"%d \n",mcells[i]);
00224         }
00225     //hanging nodes on edges
00226     for(i=0;i<nedges;i++) fprintf(stream,"%d %d %d \n",edges[2*i],edges[2*i+1],hnodes[i]);
00227 } else {//restoring original grid connectivity when me>1
00228         int lo, Ncells;
00229         double d1;
00230         stream1=fopen(grid_old,"r");
00231         scanned = fscanf(stream1,"%d \n%d \n%d \n%d \n",&lo,&i,&Ncells,&nedges);
00232         libmesh_assert_not_equal_to (scanned, EOF);
00233         fprintf(stream,"%d \n%d \n%d \n%d \n",n,N,Ncells,nedges);
00234 
00235         for(i=0;i<N;i++){//node coordinates
00236                 for(unsigned int j=0;j<n;j++) {fprintf(stream,"%e ",R[i][j]);
00237                 scanned = fscanf(stream1,"%le ",&d1);
00238                 libmesh_assert_not_equal_to (scanned, EOF);
00239         }
00240         fprintf(stream,"%d \n",mask[i]);
00241         scanned = fscanf(stream1,"%d \n",&lo);
00242         libmesh_assert_not_equal_to (scanned, EOF);
00243     }
00244         for(i=0;i<Ncells;i++){
00245                 for(unsigned int j=0;j<=3*n+n%2;j++){
00246                         scanned = fscanf(stream1,"%d ",&lo);
00247                         libmesh_assert_not_equal_to (scanned, EOF);
00248                         fprintf(stream,"%d ",lo);
00249                 }
00250                 scanned = fscanf(stream1,"\n");
00251                 libmesh_assert_not_equal_to (scanned, EOF);
00252                 fprintf(stream,"\n");
00253         }
00254         for(i=0;i<nedges;i++){
00255                 for(unsigned int j=0;j<3;j++){
00256                 scanned = fscanf(stream1,"%d ",&lo);
00257                 libmesh_assert_not_equal_to (scanned, EOF);
00258                 fprintf(stream,"%d ",lo);
00259                 }
00260         scanned = fscanf(stream1,"\n");
00261         libmesh_assert_not_equal_to (scanned, EOF);
00262         fprintf(stream,"\n");
00263         }
00264         fclose(stream1);
00265 }
00266 
00267 fclose(stream);
00268 */
00269   libMesh::out<<"Finished writegr"<<std::endl;
00270   return 0;
00271 }