example_CurlLSFEM.cpp File Reference

Example solution of a div-curl system on a hexahedral mesh using curl-conforming (edge) elements. More...

#include "TrilinosCouplings_config.h"
#include "TrilinosCouplings_Pamgen_Utils.hpp"
#include "Intrepid_FunctionSpaceTools.hpp"
#include "Intrepid_FieldContainer.hpp"
#include "Intrepid_CellTools.hpp"
#include "Intrepid_ArrayTools.hpp"
#include "Intrepid_HCURL_HEX_I1_FEM.hpp"
#include "Intrepid_HGRAD_HEX_C1_FEM.hpp"
#include "Intrepid_RealSpaceTools.hpp"
#include "Intrepid_DefaultCubatureFactory.hpp"
#include "Intrepid_Utils.hpp"
#include "Epetra_Time.h"
#include "Epetra_Map.h"
#include "Epetra_SerialComm.h"
#include "Epetra_Import.h"
#include "Epetra_Export.h"
#include "Epetra_FECrsMatrix.h"
#include "Epetra_FEVector.h"
#include "Epetra_LinearProblem.h"
#include "Teuchos_oblackholestream.hpp"
#include "Teuchos_RCP.hpp"
#include "Teuchos_ArrayRCP.hpp"
#include "Teuchos_BLAS.hpp"
#include "Teuchos_GlobalMPISession.hpp"
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_XMLParameterListHelpers.hpp"
#include "Shards_CellTopology.hpp"
#include "EpetraExt_RowMatrixOut.h"
#include "EpetraExt_MultiVectorOut.h"
#include "EpetraExt_VectorOut.h"
#include "AztecOO.h"
#include "ml_MultiLevelPreconditioner.h"
#include "ml_RefMaxwell_11_Operator.h"
#include "ml_RefMaxwell.h"
#include "ml_EdgeMatrixFreePreconditioner.h"
#include "ml_epetra_utils.h"
#include "create_inline_mesh.h"
#include "pamgen_im_exodusII_l.h"
#include "pamgen_im_ne_nemesisI_l.h"
#include "pamgen_extras.h"

Include dependency graph for example_CurlLSFEM.cpp:

Classes
struct	fecomp

Macros
#define	ABS(x)   ((x)>0?(x):-(x))

Typedefs
typedef Intrepid::FunctionSpaceTools	IntrepidFSTools
typedef Intrepid::RealSpaceTools< double >	IntrepidRSTools
typedef Intrepid::CellTools< double >	IntrepidCTools

Functions
void	TestMultiLevelPreconditioner_CurlLSFEM (char ProblemType[], Teuchos::ParameterList &MLList, Epetra_CrsMatrix &CurlCurl, Epetra_CrsMatrix &D0clean, Epetra_CrsMatrix &M0inv, Epetra_CrsMatrix &M1, Epetra_MultiVector &xh, Epetra_MultiVector &b, double &TotalErrorResidual, double &TotalErrorExactSol)
	ML Preconditioner. More...
int	evalu (double &uExact0, double &uExact1, double &uExact2, double &x, double &y, double &z)
	Exact solution evaluation. More...
double	evalDivu (double &x, double &y, double &z, double &mu)
	Divergence of exact solution. More...
int	evalCurlu (double &curlu0, double &curlu1, double &curlu2, double &x, double &y, double &z, double &mu)
	Curl of exact solution. More...
int	evalGradDivu (double &gradDivu0, double &gradDivu1, double &gradDivu2, double &x, double &y, double &z, double &mu)
	Gradient of Divergence of exact solution. More...
int	main (int argc, char *argv[])
int	Multiply_Ones (const Epetra_CrsMatrix &A, const Epetra_Vector &x, Epetra_Vector &y)
	Multiplies Ax = y, where all non-zero entries of A are replaced with the value 1.0. More...
void	solution_test (string msg, const Epetra_Operator &A, const Epetra_MultiVector &lhs, const Epetra_MultiVector &rhs, const Epetra_MultiVector &xexact, Epetra_Time &Time, double &TotalErrorExactSol, double &TotalErrorResidual)
	Compute ML solution residual. More...

Detailed Description

Example solution of a div-curl system on a hexahedral mesh using curl-conforming (edge) elements.

   This example uses the following Trilinos packages:

• Pamgen to generate a Hexahedral mesh.
• Intrepid to build the discretization matrices and right-hand side.
• Epetra to handle the global matrix and vector.
• ML to solve the linear system.

   For more details on the formulation see

• P. Bochev, K. Peterson, C. Siefert, "Analysis and Computation of Compatible Least-Squares Methods for Div-Curl Equations", SIAM J. Numerical Analysis, vol 49, pp 159-191, 2011.

        Div-Curl System:

                   curl u = g  in Omega
                    div u = h  in Omega
                    u x n = 0  on Gamma

        Corresponding discrete linear system for edge element coeficients (x):

                  (Kc + Mc*Dg*MgInv*Dg'*Mc)x = b

                  Kc    - Hcurl stiffness matrix
                  Mc    - Hcurl mass matrix
                  Dg    - Node to edge incidence matrix
                  MgInv - Hgrad mass matrix inverse
                  b     - right hand side vector

Author

Created by P. Bochev, D. Ridzal, K. Peterson, D. Hensinger, C. Siefert.

Remarks

Usage

 ./TrilinosCouplings_examples_scaling_Example_CurlLSFEM.exe  inputfile.xml


    inputfile.xml (optional)  -  xml input file containing Pamgen mesh description
                                 and material parameters for each Pamgen block,
                                 if not present code attempts to read CurlLSFEMin.xml.

Input files available in Trilinos for use with the CurlLSFEM driver:

• CurlLSFEMin.xml - basic input file with box mesh and one mesh block
• Ninja.xml - input file with distorted mesh (shaped like Ninja star) and one mesh block
• CurlLSFEMblock_in_block.xml - input file with box mesh with a center block with different material values.

Function Documentation

evalCurlu()

int evalCurlu	(	double &	curlu0,
		double &	curlu1,
		double &	curlu2,
		double &	x,
		double &	y,
		double &	z,
		double &	mu
	)

Curl of exact solution.

Parameters

curlu0	[out] first component of curl of exact solution
curlu1	[out] second component of curl of exact solution
curlu2	[out] third component of curl of exact solution
x	[in] x coordinate
y	[in] y coordinate
z	[in] z coordinate
mu	[in] material parameter

evalDivu()

double evalDivu	(	double &	x,
		double &	y,
		double &	z,
		double &	mu
	)

Divergence of exact solution.

Parameters

x	[in] x coordinate
y	[in] y coordinate
z	[in] z coordinate
mu	[in] material parameter

Returns

Value of the divergence of exact solution at (x,y,z)

evalGradDivu()

int evalGradDivu	(	double &	gradDivu0,
		double &	gradDivu1,
		double &	gradDivu2,
		double &	x,
		double &	y,
		double &	z,
		double &	mu
	)

Gradient of Divergence of exact solution.

Parameters

gradDivu0	[out] first (x) component of grad div of exact solution
gradDivu1	[out] second (y) component of grad div of exact solution
gradDivu2	[out] third (z) component of grad div of exact solution
x	[in] x coordinate
y	[in] y coordinate
z	[in] z coordinate
mu	[in] material parameter

evalu()

int evalu	(	double &	uExact0,
		double &	uExact1,
		double &	uExact2,
		double &	x,
		double &	y,
		double &	z
	)

Exact solution evaluation.

Parameters

uExact0	[out] first component of exact solution at (x,y,z)
uExact1	[out] second component of exact solution at (x,y,z)
uExact2	[out] third component of exact solution at (x,y,z)
x	[in] x coordinate
y	[in] y coordinate
z	[in] z coordinate

Multiply_Ones()

int Multiply_Ones	(	const Epetra_CrsMatrix &	A,
		const Epetra_Vector &	x,
		Epetra_Vector &	y
	)

Multiplies Ax = y, where all non-zero entries of A are replaced with the value 1.0.

Parameters

A	[in] matrix
x	[in] vector
y	[in] vector

Referenced by TestMultiLevelPreconditioner_CurlLSFEM().

solution_test()

void solution_test	(	string	msg,
		const Epetra_Operator &	A,
		const Epetra_MultiVector &	lhs,
		const Epetra_MultiVector &	rhs,
		const Epetra_MultiVector &	xexact,
		Epetra_Time &	Time,
		double &	TotalErrorExactSol,
		double &	TotalErrorResidual
	)

Compute ML solution residual.

Parameters

A	[in] discrete operator
lhs	[in] solution vector
rhs	[in] right hand side vector
Time	[in] elapsed time for output
TotalErrorResidual	[out] error residual
TotalErrorExactSol	[out] error in xh (not an appropriate measure for H(curl) basis functions)

TestMultiLevelPreconditioner_CurlLSFEM()

void TestMultiLevelPreconditioner_CurlLSFEM	(	char	ProblemType[],
		Teuchos::ParameterList &	MLList,
		Epetra_CrsMatrix &	CurlCurl,
		Epetra_CrsMatrix &	D0clean,
		Epetra_CrsMatrix &	M0inv,
		Epetra_CrsMatrix &	M1,
		Epetra_MultiVector &	xh,
		Epetra_MultiVector &	b,
		double &	TotalErrorResidual,
		double &	TotalErrorExactSol
	)

ML Preconditioner.

Parameters

ProblemType	[in] problem type
MLList	[in] ML parameter list
CurlCurl	[in] H(curl) stiffness matrix
D0clean	[in] Edge to node stiffness matrix
M0inv	[in] H(grad) mass matrix inverse
M1	[in] H(curl) mass matrix
xh	[out] solution vector
b	[in] right-hand-side vector
TotalErrorResidual	[out] error residual
TotalErrorExactSol	[out] error in xh

References Multiply_Ones().