Created
April 1, 2015 22:59
-
-
Save jwpeterson/a085afd592feca73f81f to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
// This example tests out the calculation of second derivatives | |
// Run me with: | |
// ./second_derivs-opt -d 2 -n 2 -o FIRST -ksp_monitor -pc_type hypre -pc_hypre_type boomeramg > second_derivs.log | |
// C++ include files that we need | |
#include <iostream> | |
#include <algorithm> | |
#include <math.h> | |
#include <set> | |
// Basic include file needed for the mesh functionality. | |
#include "libmesh/libmesh.h" | |
#include "libmesh/mesh.h" | |
#include "libmesh/mesh_generation.h" | |
#include "libmesh/exodusII_io.h" | |
#include "libmesh/gnuplot_io.h" | |
#include "libmesh/linear_implicit_system.h" | |
#include "libmesh/equation_systems.h" | |
#include "libmesh/fe.h" | |
#include "libmesh/quadrature_gauss.h" | |
#include "libmesh/dof_map.h" | |
#include "libmesh/sparse_matrix.h" | |
#include "libmesh/numeric_vector.h" | |
#include "libmesh/dense_matrix.h" | |
#include "libmesh/dense_vector.h" | |
#include "libmesh/perf_log.h" | |
#include "libmesh/elem.h" | |
#include "libmesh/dirichlet_boundaries.h" | |
#include "libmesh/analytic_function.h" | |
#include "libmesh/string_to_enum.h" | |
#include "libmesh/getpot.h" | |
#include "libmesh/fe_interface.h" | |
// Bring in everything from the libMesh namespace | |
using namespace libMesh; | |
// Function prototype. This is the function that will assemble | |
// the linear system for our Poisson problem. Note that the | |
// function will take the \p EquationSystems object and the | |
// name of the system we are assembling as input. From the | |
// \p EquationSystems object we have acess to the \p Mesh and | |
// other objects we might need. | |
void assemble_poisson(EquationSystems& es, | |
const std::string& system_name); | |
// Exact solution function prototype. | |
Real exact_solution (const Real x, | |
const Real y, | |
const Real z = 0.); | |
// Define a wrapper for exact_solution that will be needed below | |
void exact_solution_wrapper (DenseVector<Number>& output, | |
const Point& p, | |
const Real) | |
{ | |
output(0) = exact_solution(p(0), | |
(LIBMESH_DIM>1)?p(1):0, | |
(LIBMESH_DIM>2)?p(2):0); | |
} | |
// Begin the main program. | |
int main (int argc, char** argv) | |
{ | |
// Initialize libMesh and any dependent libaries, like in example 2. | |
LibMeshInit init (argc, argv); | |
// Declare a performance log for the main program | |
// PerfLog perf_main("Main Program"); | |
// Create a GetPot object to parse the command line | |
GetPot command_line (argc, argv); | |
// Check for proper calling arguments. | |
if (argc < 3) | |
{ | |
if (init.comm().rank() == 0) | |
libmesh_error_msg("Usage:\n" << "\t " << argv[0] << " -d 2(3)" << " -n 15"); | |
} | |
// Brief message to the user regarding the program name | |
// and command line arguments. | |
else | |
{ | |
std::cout << "Running " << argv[0]; | |
for (int i=1; i<argc; i++) | |
std::cout << " " << argv[i]; | |
std::cout << std::endl << std::endl; | |
} | |
// Read problem dimension from command line. Use int | |
// instead of unsigned since the GetPot overload is ambiguous | |
// otherwise. | |
int dim = 2; | |
if ( command_line.search(1, "-d") ) | |
dim = command_line.next(dim); | |
// Skip higher-dimensional examples on a lower-dimensional libMesh build | |
libmesh_example_requires(dim <= LIBMESH_DIM, "2D/3D support"); | |
// Create a mesh with user-defined dimension. | |
// Read number of elements from command line | |
int ps = 15; | |
if ( command_line.search(1, "-n") ) | |
ps = command_line.next(ps); | |
// Read FE order from command line | |
std::string order = "SECOND"; | |
if ( command_line.search(2, "--order", "-o") ) | |
order = command_line.next(order); | |
// An enumeration of the order string | |
Order order_enumeration = Utility::string_to_enum<Order>(order); | |
// Read FE Family from command line | |
std::string family = "LAGRANGE"; | |
if ( command_line.search(2, "--fe-family", "-f") ) | |
family = command_line.next(family); | |
// Cannot use discontinuous basis. | |
if ((family == "MONOMIAL") || (family == "XYZ")) | |
if (init.comm().rank() == 0) | |
libmesh_error_msg("ex4 currently requires a C^0 (or higher) FE basis."); | |
// An enumeration of the family string | |
FEFamily family_enumeration = Utility::string_to_enum<FEFamily>(family); | |
// Pick suitable default element type based on dim and order | |
std::string elem_type; | |
switch (order_enumeration) | |
{ | |
case FIRST: | |
{ | |
switch (dim) | |
{ | |
case 1: | |
elem_type = "EDGE2"; | |
break; | |
case 2: | |
elem_type = "QUAD4"; | |
break; | |
case 3: | |
elem_type = "HEX8"; | |
break; | |
default: | |
libmesh_error_msg("Invalid dimension: " << dim); | |
} | |
break; | |
} | |
case SECOND: | |
{ | |
switch (dim) | |
{ | |
case 1: | |
elem_type = "EDGE3"; | |
break; | |
case 2: | |
elem_type = "QUAD9"; | |
break; | |
case 3: | |
elem_type = "HEX27"; | |
break; | |
default: | |
libmesh_error_msg("Invalid dimension: " << dim); | |
} | |
break; | |
} | |
default: | |
libmesh_error_msg("Invalid order_enumeration = " << order_enumeration); | |
} | |
// Read geometric element type from command line | |
if ( command_line.search(2, "--elem-type", "-e") ) | |
elem_type = command_line.next(elem_type); | |
// An enumeration of the elem_type string | |
ElemType elem_type_enumeration = Utility::string_to_enum<ElemType>(elem_type); | |
// Should we do our finite difference check in the x, y, or z direction | |
int perturbation_index = 0; | |
if ( command_line.search(2, "--perturbation-index", "-p") ) | |
perturbation_index = command_line.next(perturbation_index); | |
if ((perturbation_index < 0) || (perturbation_index > 2)) | |
libmesh_error_msg("Error! Perturbation index must be = 0, 1, 2."); | |
// Get user's finite-differencing parameter. We do second-order central differencing, | |
// so expect the error to be O(fd_eps**2). | |
Real fd_eps = 1.e-3; | |
if ( command_line.search(1, "--fd-eps") ) | |
fd_eps = command_line.next(fd_eps); | |
// Create a mesh, with dimension to be overridden later, distributed | |
// across the default MPI communicator. | |
Mesh mesh(init.comm()); | |
Real halfwidth = dim > 1 ? 1. : 0.; | |
Real halfheight = dim > 2 ? 1. : 0.; | |
// Use the MeshTools::Generation mesh generator to create a uniform | |
// grid on the square [-1,1]^D. For non-EDGE/QUAD/HEX | |
// elem_type_enumerations, build_cube() will produce meshes with more | |
// than one element -- e.g. 6 for pyramids and 24 for tets. | |
MeshTools::Generation::build_cube (mesh, | |
ps, | |
(dim>1) ? ps : 0, | |
(dim>2) ? ps : 0, | |
-1., 1., | |
-halfwidth, halfwidth, | |
-halfheight, halfheight, | |
elem_type_enumeration); | |
// Print information about the mesh to the screen. | |
mesh.print_info(); | |
// Node 0 will always be at (-1,-1) in 2D and (-1,-1,-1) in 3D. It | |
// should be fairly safe and general to displace this node in order | |
// to verify that the second derivatives are calculated correctly on | |
// non-affine elements. | |
// libMesh::out << mesh.point(0) << std::endl; | |
Node* node = mesh.node_ptr(0); | |
if (dim==2) | |
{ | |
(*node)(0) = -1.1; | |
(*node)(1) = -1.1; | |
} | |
else if (dim==3) | |
{ | |
(*node)(0) = -1.1; | |
(*node)(1) = -1.1; | |
(*node)(2) = -1.1; | |
} | |
// Create an equation systems object. | |
EquationSystems equation_systems (mesh); | |
// Store user's requested perturbation_index in the ES object. | |
equation_systems.parameters.set<int>("perturbation_index") = perturbation_index; | |
// Set the user's requested finite-difference parameter | |
equation_systems.parameters.set<Real>("fd_eps") = fd_eps; | |
// Declare the system and its variables. | |
// Create a system named "Poisson" | |
LinearImplicitSystem& system = | |
equation_systems.add_system<LinearImplicitSystem> ("Poisson"); | |
// Add the variable "u" to "Poisson". "u" | |
// will be approximated using second-order approximation. | |
unsigned int u_var = system.add_variable("u", | |
order_enumeration, | |
family_enumeration); | |
// Give the system a pointer to the matrix assembly | |
// function. | |
system.attach_assemble_function (assemble_poisson); | |
// Construct a Dirichlet boundary condition object | |
// Indicate which boundary IDs we impose the BC on | |
// We either build a line, a square or a cube, and | |
// here we indicate the boundaries IDs in each case | |
std::set<boundary_id_type> boundary_ids; | |
// the dim==1 mesh has two boundaries with IDs 0 and 1 | |
boundary_ids.insert(0); | |
boundary_ids.insert(1); | |
// the dim==2 mesh has four boundaries with IDs 0, 1, 2 and 3 | |
if(dim>=2) | |
{ | |
boundary_ids.insert(2); | |
boundary_ids.insert(3); | |
} | |
// the dim==3 mesh has four boundaries with IDs 0, 1, 2, 3, 4 and 5 | |
if(dim==3) | |
{ | |
boundary_ids.insert(4); | |
boundary_ids.insert(5); | |
} | |
// Create a vector storing the variable numbers which the BC applies to | |
std::vector<unsigned int> variables(1); | |
variables[0] = u_var; | |
// Create an AnalyticFunction object that we use to project the BC | |
// This function just calls the function exact_solution via exact_solution_wrapper | |
AnalyticFunction<> exact_solution_object(exact_solution_wrapper); | |
DirichletBoundary dirichlet_bc(boundary_ids, | |
variables, | |
&exact_solution_object); | |
// We must add the Dirichlet boundary condition _before_ | |
// we call equation_systems.init() | |
system.get_dof_map().add_dirichlet_boundary(dirichlet_bc); | |
// Initialize the data structures for the equation system. | |
equation_systems.init(); | |
// Print information about the system to the screen. | |
equation_systems.print_info(); | |
mesh.print_info(); | |
// Solve the system "Poisson", just like example 2. | |
system.solve(); | |
// After solving the system write the solution | |
// to a GMV-formatted plot file. | |
if(dim == 1) | |
{ | |
GnuPlotIO plot(mesh,"Introduction Example 4, 1D",GnuPlotIO::GRID_ON); | |
plot.write_equation_systems("gnuplot_script",equation_systems); | |
} | |
#ifdef LIBMESH_HAVE_EXODUS_API | |
else | |
{ | |
ExodusII_IO (mesh).write_equation_systems ((dim == 3) ? | |
"out_3.e" : "out_2.e",equation_systems); | |
} | |
#endif // #ifdef LIBMESH_HAVE_EXODUS_API | |
// All done. | |
return 0; | |
} | |
// We now define the matrix assembly function for the | |
// Poisson system. We need to first compute element | |
// matrices and right-hand sides, and then take into | |
// account the boundary conditions. | |
void assemble_poisson(EquationSystems& es, | |
const std::string& system_name) | |
{ | |
// It is a good idea to make sure we are assembling | |
// the proper system. | |
libmesh_assert_equal_to (system_name, "Poisson"); | |
// Declare a performance log. Give it a descriptive | |
// string to identify what part of the code we are | |
// logging, since there may be many PerfLogs in an | |
// application. | |
PerfLog perf_log ("Matrix Assembly"); | |
// Get a constant reference to the mesh object. | |
const MeshBase& mesh = es.get_mesh(); | |
// The dimension that we are running | |
const unsigned int dim = mesh.mesh_dimension(); | |
// Get a reference to the LinearImplicitSystem we are solving | |
LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("Poisson"); | |
// A reference to the \p DofMap object for this system. The \p DofMap | |
// object handles the index translation from node and element numbers | |
// to degree of freedom numbers. We will talk more about the \p DofMap | |
// in future examples. | |
const DofMap& dof_map = system.get_dof_map(); | |
// Get a constant reference to the Finite Element type | |
// for the first (and only) variable in the system. | |
FEType fe_type = dof_map.variable_type(0); | |
// Build a Finite Element object of the specified type. Since the | |
// \p FEBase::build() member dynamically creates memory we will | |
// store the object as an \p UniquePtr<FEBase>. This can be thought | |
// of as a pointer that will clean up after itself. | |
UniquePtr<FEBase> fe (FEBase::build(dim, fe_type)); | |
// For QUAD/HEX | |
// THIRD=2x2x...x2 | |
// FIFTH=3x3x...x3 | |
Order qorder = THIRD; | |
// A 5th order Gauss quadrature rule for numerical integration. | |
QGauss qrule (dim, qorder); | |
// Tell the finite element object to use our quadrature rule. | |
fe->attach_quadrature_rule (&qrule); | |
// Declare a special finite element object for | |
// boundary integration. | |
UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type)); | |
// Boundary integration requires one quadraure rule, | |
// with dimensionality one less than the dimensionality | |
// of the element. | |
QGauss qface(dim-1, qorder); | |
// Tell the finte element object to use our | |
// quadrature rule. | |
fe_face->attach_quadrature_rule (&qface); | |
// Here we define some references to cell-specific data that | |
// will be used to assemble the linear system. | |
// We begin with the element Jacobian * quadrature weight at each | |
// integration point. | |
const std::vector<Real>& JxW = fe->get_JxW(); | |
// The physical XY locations of the quadrature points on the element. | |
// These might be useful for evaluating spatially varying material | |
// properties at the quadrature points. | |
const std::vector<Point>& q_point = fe->get_xyz(); | |
// The element shape functions evaluated at the quadrature points. | |
const std::vector<std::vector<Real> >& phi = fe->get_phi(); | |
// The element shape function gradients evaluated at the quadrature | |
// points. | |
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi(); | |
// The element shape function second derivatives evaluated at the | |
// quadrature points. | |
const std::vector<std::vector<RealTensor> >& d2phi = fe->get_d2phi(); | |
// Define data structures to contain the element matrix | |
// and right-hand-side vector contribution. Following | |
// basic finite element terminology we will denote these | |
// "Ke" and "Fe". More detail is in example 3. | |
DenseMatrix<Number> Ke; | |
DenseVector<Number> Fe; | |
// This vector will hold the degree of freedom indices for | |
// the element. These define where in the global system | |
// the element degrees of freedom get mapped. | |
std::vector<dof_id_type> dof_indices; | |
// Now we will loop over all the elements in the mesh. | |
// We will compute the element matrix and right-hand-side | |
// contribution. See example 3 for a discussion of the | |
// element iterators. | |
MeshBase::const_element_iterator el = mesh.active_local_elements_begin(); | |
const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end(); | |
for ( ; el != end_el; ++el) | |
{ | |
// Start logging the shape function initialization. | |
// This is done through a simple function call with | |
// the name of the event to log. | |
perf_log.push("elem init"); | |
// Store a pointer to the element we are currently | |
// working on. This allows for nicer syntax later. | |
const Elem* elem = *el; | |
std::cout << "\nAssembling element: " << elem->id() << std::endl; | |
// Get the degree of freedom indices for the | |
// current element. These define where in the global | |
// matrix and right-hand-side this element will | |
// contribute to. | |
dof_map.dof_indices (elem, dof_indices); | |
// Compute the element-specific data for the current | |
// element. This involves computing the location of the | |
// quadrature points (q_point) and the shape functions | |
// (phi, dphi) for the current element. | |
fe->reinit (elem); | |
// Zero the element matrix and right-hand side before | |
// summing them. We use the resize member here because | |
// the number of degrees of freedom might have changed from | |
// the last element. Note that this will be the case if the | |
// element type is different (i.e. the last element was a | |
// triangle, now we are on a quadrilateral). | |
Ke.resize (dof_indices.size(), | |
dof_indices.size()); | |
Fe.resize (dof_indices.size()); | |
// Stop logging the shape function initialization. | |
// If you forget to stop logging an event the PerfLog | |
// object will probably catch the error and abort. | |
perf_log.pop("elem init"); | |
// Now we will build the element matrix. This involves | |
// a double loop to integrate the test funcions (i) against | |
// the trial functions (j). | |
// | |
// We have split the numeric integration into two loops | |
// so that we can log the matrix and right-hand-side | |
// computation seperately. | |
// | |
// Now start logging the element matrix computation | |
perf_log.push ("Ke"); | |
for (unsigned int qp=0; qp<qrule.n_points(); qp++) | |
for (unsigned int i=0; i<phi.size(); i++) | |
for (unsigned int j=0; j<phi.size(); j++) | |
Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]); | |
// Stop logging the matrix computation | |
perf_log.pop ("Ke"); | |
// Now we build the element right-hand-side contribution. | |
// This involves a single loop in which we integrate the | |
// "forcing function" in the PDE against the test functions. | |
// | |
// Start logging the right-hand-side computation | |
perf_log.push ("Fe"); | |
for (unsigned int qp=0; qp<qrule.n_points(); qp++) | |
{ | |
// fxy is the forcing function for the Poisson equation. | |
// In this case we set fxy to be a finite difference | |
// Laplacian approximation to the (known) exact solution. | |
// | |
// We will use the second-order accurate FD Laplacian | |
// approximation, which in 2D on a structured grid is | |
// | |
// u_xx + u_yy = (u(i-1,j) + u(i+1,j) + | |
// u(i,j-1) + u(i,j+1) + | |
// -4*u(i,j))/h^2 | |
// | |
// Since the value of the forcing function depends only | |
// on the location of the quadrature point (q_point[qp]) | |
// we will compute it here, outside of the i-loop | |
const Real x = q_point[qp](0); | |
#if LIBMESH_DIM > 1 | |
const Real y = q_point[qp](1); | |
#else | |
const Real y = 0.; | |
#endif | |
#if LIBMESH_DIM > 2 | |
const Real z = q_point[qp](2); | |
#else | |
const Real z = 0.; | |
#endif | |
const Real eps = 1.e-3; | |
const Real uxx = (exact_solution(x-eps,y,z) + | |
exact_solution(x+eps,y,z) + | |
-2.*exact_solution(x,y,z))/eps/eps; | |
const Real uyy = (exact_solution(x,y-eps,z) + | |
exact_solution(x,y+eps,z) + | |
-2.*exact_solution(x,y,z))/eps/eps; | |
const Real uzz = (exact_solution(x,y,z-eps) + | |
exact_solution(x,y,z+eps) + | |
-2.*exact_solution(x,y,z))/eps/eps; | |
Real fxy = - (uxx + ((dim==1) ? 0. : uyy) + ((dim==2) ? 0. : uzz)); | |
// Add the RHS contribution | |
for (unsigned int i=0; i<phi.size(); i++) | |
Fe(i) += JxW[qp]*fxy*phi[i][qp]; | |
// libMesh::out << std::endl; | |
// // Real U_xx=0., U_xy=0.; | |
// for (unsigned int i=0; i<d2phi.size(); i++) | |
// { | |
// // U_st = sum u_i phi_{i,st} | |
// // U_xx += system.current_solution(dof_indices[i])*d2phi[i][qp](0,0); | |
// // U_xy += system.current_solution(dof_indices[i])*d2phi[i][qp](0,1); | |
// | |
// // At each qp, the second derivative is a tensor of values | |
// // xx xy xz | |
// // yx yy yz | |
// // zx zy zz | |
// if (i==0) | |
// std::cout << "x_qp = " << q_point[qp] << std::endl; | |
// std::cout << "D^2 phi[i=" << i << "]= " << std::endl; | |
// std::cout << d2phi[i][qp]; | |
// } | |
// libMesh::out << "U_xx[" << qp << "]= " << U_xx << std::endl; | |
// libMesh::out << "U_xy[" << qp << "]= " << U_xy << std::endl; | |
} | |
// Stop logging the right-hand-side computation | |
perf_log.pop ("Fe"); | |
// Make a copy of shape function quantities before we reinitialize the quadrature rule at new points. | |
const std::vector<std::vector<Real> > center_phi = phi; | |
const std::vector<std::vector<RealGradient> > center_dphi = dphi; | |
const std::vector<std::vector<RealTensor> > center_d2phi = d2phi; | |
// 0=perturb in x direction, 1=perturb in y, etc | |
const int perturbation_index = es.parameters.get<int>("perturbation_index"); | |
// Finite differencing parameter, as requested by user | |
const Real fd_eps = es.parameters.get<Real>("fd_eps"); | |
{ | |
// Forward perturbation vector | |
Point delta; | |
delta(perturbation_index) = fd_eps; | |
// Perturb quadrature points in *physical* space, inverse_map() | |
// them back to reference space, in preparation for doing a | |
// finite difference calculation of various element data to | |
// check second derivative quantities. | |
std::vector<Point> forward_perturbed_physical_points(qrule.n_points()); | |
for (unsigned int qp=0; qp<qrule.n_points(); qp++) | |
{ | |
// Perturb by fd_eps | |
forward_perturbed_physical_points[qp] = q_point[qp] + delta; | |
} | |
// Inverse map the perturbed point back to reference space | |
std::vector<Point> forward_perturbed_reference_points(qrule.n_points()); | |
FEInterface::inverse_map(dim, | |
fe_type, | |
elem, | |
forward_perturbed_physical_points, | |
forward_perturbed_reference_points); | |
// Print before/after | |
// std::cout << "\nOriginal quadrature points." << std::endl; | |
// for (unsigned int qp=0; qp<qrule.n_points(); qp++) | |
// std::cout << qrule.get_points()[qp] << std::endl; | |
// | |
// std::cout << "\nForward perturbed quadrature points." << std::endl; | |
// for (unsigned int qp=0; qp<qrule.n_points(); qp++) | |
// std::cout << forward_perturbed_reference_points[qp] << std::endl; | |
// Reinitialize the FE at the perturbed quadrature points | |
fe->reinit (elem, &forward_perturbed_reference_points); | |
} | |
// Save the forward-perturbed shape function values | |
const std::vector<std::vector<Real> > forward_phi = phi; | |
const std::vector<std::vector<RealGradient> > forward_dphi = dphi; | |
const std::vector<std::vector<RealTensor> > forward_d2phi = d2phi; | |
{ | |
// Reverse perturbation vectors. Note: because I already perturbed | |
// forward by +fd_eps, I have to perturb backward by 2 fd_eps to get | |
// to a backward perturbed state.. | |
Point delta; | |
delta(perturbation_index) = -2.*fd_eps; | |
// Perturb quadrature points in *physical* space, inverse_map() | |
// them back to reference space, in preparation for doing a | |
// finite difference calculation of various element data to | |
// check second derivative quantities. | |
std::vector<Point> reverse_perturbed_physical_points(qrule.n_points()); | |
for (unsigned int qp=0; qp<qrule.n_points(); qp++) | |
{ | |
// Perturb by fd_eps | |
reverse_perturbed_physical_points[qp] = q_point[qp] + delta; | |
} | |
// Inverse map the perturbed point back to reference space | |
std::vector<Point> reverse_perturbed_reference_points(qrule.n_points()); | |
FEInterface::inverse_map(dim, | |
fe_type, | |
elem, | |
reverse_perturbed_physical_points, | |
reverse_perturbed_reference_points); | |
// Reinitialize the FE at the perturbed quadrature points | |
fe->reinit (elem, &reverse_perturbed_reference_points); | |
} | |
// Save the reverse-perturbed shape function values | |
const std::vector<std::vector<Real> > reverse_phi = phi; | |
const std::vector<std::vector<RealGradient> > reverse_dphi = dphi; | |
const std::vector<std::vector<RealTensor> > reverse_d2phi = d2phi; | |
// Look at values of various things at perturbed quadrature points | |
for (unsigned int qp=0; qp<qrule.n_points(); qp++) | |
{ | |
for (unsigned int i=0; i<dphi.size(); i++) | |
{ | |
// First derivative quantities | |
Real abs_err_1st[3] = {0., 0., 0.}; // d/dx=0, d/dy=1, d/dz=2 | |
Real fd_approx_1st[3] = {0., 0., 0.}; // d/dx=0, d/dy=1, d/dz=2 | |
Real fe_1st[3] = {0., 0., 0.}; // d/dx=0, d/dy=1, d/dz=2 | |
const char* deriv_strings_1st[3] = {"x", "y", "z"}; | |
// Second derivative quantities | |
Real abs_err[9] = | |
{ | |
0., 0., 0., // xx=0, xy=1, xz=2 | |
0., 0., 0., // yx=3, yy=4, yz=5 | |
0., 0., 0. // zx=6, zy=7, zz=8 | |
}; | |
Real fd_approx_2nd[9] = | |
{ | |
0., 0., 0., // xx=0, xy=1, xz=2 | |
0., 0., 0., // yx=3, yy=4, yz=5 | |
0., 0., 0. // zx=6, zy=7, zz=8 | |
}; | |
// The second derivative, as computed by the FE object. | |
Real fe_2nd[9] = | |
{ | |
0., 0., 0., // xx=0, xy=1, xz=2 | |
0., 0., 0., // yx=3, yy=4, yz=5 | |
0., 0., 0. // zx=6, zy=7, zz=8 | |
}; | |
const char* deriv_strings[9] = {"xx", "xy", "xz", "yx", "yy", "yz", "zx", "zy", "zz"}; | |
// (Depending on perturbation_index!) compute: | |
// u_x \approx (u(x+fd_eps) - u(x-fd_eps))/2/fd_eps | |
// u_y \approx (u(y+fd_eps) - u(y-fd_eps))/2/fd_eps | |
// u_z \approx (u(z+fd_eps) - u(z-fd_eps))/2/fd_eps | |
fd_approx_1st[perturbation_index] = (forward_phi[i][qp] - reverse_phi[i][qp])/2./fd_eps; | |
fe_1st[perturbation_index] = center_dphi[i][qp](perturbation_index); | |
abs_err_1st[perturbation_index] = std::abs(fd_approx_1st[perturbation_index] - fe_1st[perturbation_index]); | |
// Compute whichever 3 of the following second derivatives makes sense for the current perturbation_index: | |
// u_{xy} = u_{yx} \approx (u_x(y+fd_eps) - u_x(y-fd_eps))/2/fd_eps | |
// u_{yy} \approx (u_y(y+fd_eps) - u_y(y-fd_eps))/2/fd_eps | |
// u_{zy} = u_{yz} \approx (u_z(y+fd_eps) - u_z(y-fd_eps))/2/fd_eps | |
// | |
// u_{xy} = u_{yx} \approx (u_x(y+fd_eps) - u_x(y-fd_eps))/2/fd_eps | |
// u_{yy} \approx (u_y(y+fd_eps) - u_y(y-fd_eps))/2/fd_eps | |
// u_{zy} = u_{yz} \approx (u_z(y+fd_eps) - u_z(y-fd_eps))/2/fd_eps | |
// | |
// u_{xz} = u_{zx} \approx (u_x(z+fd_eps) - u_x(z-fd_eps))/2/fd_eps | |
// u_{zy} = u_{yz} \approx (u_y(z+fd_eps) - u_y(z-fd_eps))/2/fd_eps | |
// u_{zz} \approx (u_z(z+fd_eps) - u_z(z-fd_eps))/2/fd_eps | |
for (unsigned int j=0; j<3; ++j) | |
{ | |
fe_2nd[j+3*perturbation_index] = center_d2phi[i][qp](perturbation_index, j); | |
fd_approx_2nd[j+3*perturbation_index] = (forward_dphi[i][qp](j) - reverse_dphi[i][qp](j))/2./fd_eps; | |
abs_err[j+3*perturbation_index] = std::abs(fd_approx_2nd[j+3*perturbation_index] - fe_2nd[j+3*perturbation_index]); | |
} | |
// 1st derivative error messages | |
std::ios_base::fmtflags cout_flags = std::cout.flags(); | |
std::cout << std::scientific << std::setprecision(16); | |
for (unsigned j=0; j<3; ++j) | |
{ | |
// if (abs_err_1st[j] > 1.e-10) // Print if error is above tolerance | |
if (abs_err_1st[j] > 0.) // Print any non-zero error | |
{ | |
std::cout << "i=" << i << ", qp=" << qp << " = " | |
<< "FE Approx u_" << deriv_strings_1st[j] << "=" << fe_1st[j] | |
<< ", FD Approx u_" << deriv_strings_1st[j] << "=" << fd_approx_1st[j] | |
<< ", Abs. Error in u_" << deriv_strings_1st[j] << "=" << abs_err_1st[j] << std::endl; | |
} | |
} | |
// 2nd derivative error messages | |
for (unsigned j=0; j<9; ++j) | |
{ | |
// if (abs_err[j] > 1.e-10) // Print if error is above tolerance | |
if (abs_err[j] > 0.) // Print any non-zero error | |
{ | |
std::cout << "i=" << i << ", qp=" << qp << " = " | |
<< "FE Approx u_" << deriv_strings[j] << "=" << fe_2nd[j] | |
<< ", FD Approx u_" << deriv_strings[j] << "=" << fd_approx_2nd[j] | |
<< ", Abs. Error in u_" << deriv_strings[j] << "=" << abs_err[j] << std::endl; | |
// Actually *stop* executing when we find such an error? | |
// libmesh_error_msg("Error between FD and FE approximation to second derivative!"); | |
} | |
} | |
// Restore formatting flags | |
std::cout.flags(cout_flags); | |
} // end for i | |
} // end for qp | |
// If this assembly program were to be used on an adaptive mesh, | |
// we would have to apply any hanging node constraint equations | |
// Also, note that here we call heterogenously_constrain_element_matrix_and_vector | |
// to impose a inhomogeneous Dirichlet boundary conditions. | |
dof_map.heterogenously_constrain_element_matrix_and_vector (Ke, Fe, dof_indices); | |
// The element matrix and right-hand-side are now built | |
// for this element. Add them to the global matrix and | |
// right-hand-side vector. The \p SparseMatrix::add_matrix() | |
// and \p NumericVector::add_vector() members do this for us. | |
// Start logging the insertion of the local (element) | |
// matrix and vector into the global matrix and vector | |
perf_log.push ("matrix insertion"); | |
system.matrix->add_matrix (Ke, dof_indices); | |
system.rhs->add_vector (Fe, dof_indices); | |
// Start logging the insertion of the local (element) | |
// matrix and vector into the global matrix and vector | |
perf_log.pop ("matrix insertion"); | |
} | |
// That's it. We don't need to do anything else to the | |
// PerfLog. When it goes out of scope (at this function return) | |
// it will print its log to the screen. Pretty easy, huh? | |
} | |
Real exact_solution (const Real x, | |
const Real y, | |
const Real /*z*/) | |
{ | |
return x*y; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment