//LIC// ====================================================================
//LIC// This file forms part of oomph-lib, the object-oriented,
//LIC// multi-physics finite-element library, available
//LIC// at http://www.oomph-lib.org.
//LIC//
//LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
//LIC//
//LIC// This library is free software; you can redistribute it and/or
//LIC// modify it under the terms of the GNU Lesser General Public
//LIC// License as published by the Free Software Foundation; either
//LIC// version 2.1 of the License, or (at your option) any later version.
//LIC//
//LIC// This library is distributed in the hope that it will be useful,
//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
//LIC// Lesser General Public License for more details.
//LIC//
//LIC// You should have received a copy of the GNU Lesser General Public
//LIC// License along with this library; if not, write to the Free Software
//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
//LIC// 02110-1301 USA.
//LIC//
//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
//LIC//
//LIC//====================================================================
//Demonstration code to solve a one-dimensional Poisson equation
//using the OOMPH-LIB base classes only
//Include std C++ IO library
#include
//Include the OOMPH-LIB base classes
#include "generic.h"
using namespace std;
using namespace oomph;
//------------------------GEOMETRIC ELEMENT----------------------------------
//===========================================================================
/// A "Geometric" Element consists of Nodes only and contains all the
/// information necessary to find the value of one, or more, global
/// coordinates at the nodes and at any point within the element. It should
/// also specify a spatial integration scheme that can be used to
/// calculate the size of the element.
///
/// A linear, one dimensional geometric element consists of two Nodes and
/// uses linear interpolation between the Nodes. The integration scheme is
/// a two-point Gaussian quadrature.
//===========================================================================
class TwoNodeGeometricElement : public FiniteElement
{
private:
/// Integration scheme that will be used to integrate over the element.
/// Simple Gaussian quadrature in one dimension, with two Gauss points.
static Gauss<1,2> Default_spatial_integration_scheme;
public:
/// Element constructor; resize Node_pt to the number of nodes (two)
/// and set the integration scheme.
TwoNodeGeometricElement()
{
//Linear interpolation requires two Nodes per element.
//In fact, calling this function merely provides storage for
//the pointers to the Nodes and initialises the pointers to NULL.
//The Nodes themselves are created during the mesh generation
//process by the functions FiniteElement::construct_node(...)
//which stores the pointers to the newly created Nodes
//in the element's own internal storage.
this->set_n_node(2);
//The element is one-dimensional
this->set_dimension(1);
//Set the pointer to the spatial integration scheme
set_integration_scheme(&Default_spatial_integration_scheme);
}
/// Define the shape functions for interpolation across the element.
void shape(const Vector &s, Shape &psi) const
{
//There are two shape functions (one per node)
//In terms of the local coordinate, s[0]:
//Node 0 is at s[0] = -1.0
//Node 1 is at s[0] = 1.0
//psi[0] takes the value one at node 0 and zero at node 1
psi[0] = 0.5*(1.0 - s[0]);
//psi[1] takes the value one at node 1 and zero at node 0
psi[1] = 0.5*(1.0 + s[0]);
}
/// Define the derivatives of the shape functions w.r.t. the local coordinate
void dshape_local(const Vector &s, Shape &psi, DShape &dpsids) const
{
//Call the shape functions
shape(s,psi);
//The derivative of psi[0] wrt s[0] is -0.5
dpsids(0,0) = -0.5;
//The derivative of psi[0] wrt s[0] is 0.5
dpsids(1,0) = 0.5;
}
/// The element consists of two nodes along its one (and only) edge
unsigned nnode_1d() const {return 2;}
};
//-----------------------SPECIFIC FINITE ELEMENT---------------------------
//=========================================================================
/// A specific Element needs to implement the functions that compute
/// the elemental contributions to the global residual vector and the
/// global Jacobian matrix formed when using Newton's method.
/// By using calls to the interpolation scheme defined in an
/// an underlying "Geometric" Element, the equations can be written
/// generically and Elements that solve the same equations,
/// but using different orders of interpolation may be easily
/// constructed, using a different "geometric" base class.
///
/// A linear, Poisson element solves a one-dimensional Poisson equation
/// using the integration and interpolation scheme defined in the
/// above "Geometric" element.
//=========================================================================
class TwoNodePoissonElement : public TwoNodeGeometricElement
{
public:
/// Constructor: Initialise the sign of the source function to +1
TwoNodePoissonElement()
{
// Initialise sign of source function
Sign=1;
}
/// Define a specific source function. For greater generality, and
/// inclusion in a library, this could be defined as a source function
/// pointer, that would then be set externally.
double f(const double &x) const
{
return double(Sign)*30.0*sin(sqrt(30.0)*x);
}
/// Access function to the sign in the source function
int& sign() {return Sign;}
/// Define an access function to the first data value stored
/// at each node. In a more general "Equation" element, such abstraction
/// is essential, because different Elements will store the same variables
/// in different locations.
double u(const unsigned &n) {return node_pt(n)->value(0);}
/// For the Poisson equation, only one value is stored at each node
unsigned required_nvalue(const unsigned &n) const {return 1;}
/// Calculate the elemental contributions to the global
/// residual vector for the weak form of the Poisson equation
void get_residuals(Vector &residuals)
{
//Find the number of degrees of freedom (unpinned values) in the element
unsigned n_dof = ndof();
//Initialise all the residuals to zero
for(unsigned i=0;i s(1);
//Find the number of integration points in the underlying
//geometric element's integration scheme
unsigned n_intpt = integral_pt()->nweight();
//Loop over the integration points
for(unsigned ipt=0;iptknot(ipt,0);
//Find the weight of the integration scheme at this knot point
double w = integral_pt()->weight(ipt);
//Find the shape functions and their derivatives at the knot point.
//This function is implemented in FiniteElement.
//It also returns the Jacobian of the mapping from local to
//global coordinates.
double J = dshape_eulerian(s,psi,dpsidx);
//Premultiply the weight and the Jacobian
double W = w*J;
//Allocate storage for the value of the field variable u,
//its derivative and the global position at the knot point.
//Initialise them all to zero.
double interpolated_x=0.0, interpolated_u=0.0, interpolated_dudx=0.0;
//Calculate the interpolated values by looping over the shape
//functions and summing the appropriate contributions
for(unsigned n=0;n= 0)
{
//Add body force/source term here
residuals[local_eqn_number] += source*psi[l]*W;
//Add the Poisson bit itself
residuals[local_eqn_number] += interpolated_dudx*dpsidx(l,0)*W;
}
}
} //End of loop over the integration points
} //End of function
/// Calculate the elemental contribution to the global residual
/// vector and to the Jacobian matrix dR_{i}/du_{j} used in the Newton method
void get_jacobian(Vector &residuals, DenseMatrix &jacobian)
{
//First, calculate the residuals
get_residuals(residuals);
//Find the number of degrees of freedom (unpinned values) in the element
unsigned n_dof = ndof();
//Initialise all entries of the Jacobian matrix to zero
for(unsigned i=0;i s(1);
//Find the number of integration points in the underlying
//geometric element's integration scheme
unsigned n_intpt = integral_pt()->nweight();
//Loop over the integration points
for(unsigned ipt=0;iptknot(ipt,0);
//Find the weight of the integration scheme at this knot point
double w = integral_pt()->weight(ipt);
//Find the shape functions and their derivatives at the knot point.
//This function is implemented in FiniteElement.
//It also returns the Jacobian of the mapping from local to
//global coordinates.
double J = dshape_eulerian(s,psi,dpsidx);
//Premultiply the weight and the Jacobian
double W = w*J;
//ASSEMBLE THE JACOBIAN TERMS
//Loop over the test (shape) functions
for(unsigned l=0;l= 0)
{
//Loop over the degrees of freedom
for(unsigned l2=0;l2= 0)
{
//Add the contribution to the Jacobian
jacobian(local_eqn_number,local_dof_number) +=
dpsidx(l2,0)*dpsidx(l,0)*W;
}
}
}
}
} //End of loop over the integration points
} //End of function
//Define an output function for the element
void output(ostream &output)
{
//Read out the number of nodes in the element
unsigned n_node = nnode();
//Loop over the nodes and print out the global coordinate
//and value of the field variable, u, at each node
for(unsigned n=0;n TwoNodeGeometricElement::Default_spatial_integration_scheme;
//----------------------ONE DIMENSIONAL (LINE) MESH---------------------------
//============================================================================
/// A simple one dimensional mesh: uniformly spaced nodes in the domain x=0,1
//============================================================================
template
class OneDimMesh : public Mesh
{
public:
/// Mesh Constructor. The argument is the desired number of elements
OneDimMesh(const unsigned &n_element)
{
//Resize the vector of pointers to elements: there are n_element elements
Element_pt.resize(n_element);
//Construct the first element (Note the use of the template parameter)
Element_pt[0] = new ELEMENT;
//Construct the first node and add it to the Mesh::Node_pt vector
//Note: The FiniteElement::construct_boundary_node(j) function
//builds the element's j-th local node, and provides the functionality
//that allows it to be located on a Mesh boundary -- essentially this
//involves allocating additional storage to the Node.
//The function obtains the Node's
//characteristics (e.g. its spatial dimension, the number of
//values to be stored, etc) from various virtual FiniteElement
//member functions, such as FiniteElement::required_nvalue().
//FiniteElement::construct_boundary_node(...) also
//stores a pointer to the newly created Node in the element's own
//Node_pt vector.
//Finally, the function returns the pointer to the
//newly created Node, so that it can be stored in the Mesh's Node_pt
//vector, as done here:
Node_pt.push_back(finite_element_pt(0)->construct_boundary_node(0));
//Find the number of nodes per element (N.B. all elements are identical
//so we can determine this value once and for all).
unsigned n_node = finite_element_pt(0)->nnode();
//Loop over the remaning nodes of the first element
for(unsigned n=1;nconstruct_node(n));
}
//Loop over the remaining elements apart from the last
for(unsigned e=1;e<(n_element-1);e++)
{
//Construct the e-th element
Element_pt[e] = new ELEMENT;
//The first local node of the e-th element is the last local node
//of the (e-1)-th element. We MUST NOT construct the node twice.
//Instead, we set the pointer in the e-th element to point to the
//previously created node in the (e-1)-th element.
finite_element_pt(e)->node_pt(0) =
finite_element_pt(e-1)->node_pt(n_node-1);
//Loop over the remaining nodes of the e-th element
for(unsigned n=1;nconstruct_node(n));
}
} //End of loop over elements
//Construct the final element
Element_pt[n_element-1] = new ELEMENT;
//The first local node of the final element is the last local node
//of the penultimate element. We MUST NOT construct the node twice.
//Instead, we set the pointer in the final element to point to the
//previously created node in the penultimate element.
finite_element_pt(n_element-1)->node_pt(0) =
finite_element_pt(n_element-2)->node_pt(n_node-1);
//Loop over the remaining central nodes of the final element
for(unsigned n=1;n<(n_node-1);n++)
{
//Construct the next node and add it to the Mesh::Node_pt vector
//Note that these interior nodes need not (and should not)
//be boundary nodes, so they are created using the construct_node
//function()
Node_pt.push_back(finite_element_pt(n_element-1)->construct_node(n));
}
//Construct the final node and add it to the Mesh::Node_pt vector.
//This node will be located on a boundary, and hence we use
//the construct_boundary_node function.
Node_pt.push_back(finite_element_pt(n_element-1)
->construct_boundary_node(n_node-1));
//We've now created all the nodes -- let's set their positions:
//Find the total number of nodes
unsigned n_global_node = nnode();
//Loop over all nodes
for(unsigned n=0;nx(0) = double(n)/double(n_global_node-1);
}
//Set the boundary data:
//There are two boundaries in this mesh
set_nboundary(2);
//Boundary 0 contains the first node in the mesh:
add_boundary_node(0,Node_pt[0]);
//Boundary 1 contains the final node in the mesh:
add_boundary_node(1,Node_pt[n_global_node-1]);
} // End of constructor
}; // End of OneDimMesh class.
//-----------------------PROBLEM CLASS-------------------------------------
//==========================================================================
/// Define the Problem which creates the mesh, applies the
/// boundary conditions, and assigns equation numbers.
//==========================================================================
class DemoPoissonProblem : public Problem
{
public:
/// Problem constructor: Pass the sign of the source function (default
/// is +1)
DemoPoissonProblem(const int& sign=1) : Sign(sign)
{
//Create a OneDimMesh Mesh object and set it to be the problem's mesh.
//The element type, TwoNodePoissonElement, is passed as a template
//parameter to the mesh. The argument to the constructor indicates
//the number of elements in the mesh.
Problem::mesh_pt() = new OneDimMesh(10);
//Pin the unknowns at the ends of the 1D domain:
//The 1D mesh has 2 boundaries, each of which contains a single node;
//the nodes on the boundary are available from Mesh::boundary_node_pt(...)
//Pin the single nodal value at the single node on mesh
//boundary 0 (= the left domain boundary at x=0)
mesh_pt()->boundary_node_pt(0,0)->pin(0);
//Pin the single nodal value at the single node on mesh
//boundary 1 (= the right domain boundary at x=1)
mesh_pt()->boundary_node_pt(1,0)->pin(0);
// All values are initialised to zero. This is consistent with the
// boundary condition at x=0 and no further action is required
// at that node.
// Apply the boundary condition at x=1: u(x=1)=-/+1
mesh_pt()->boundary_node_pt(1,0)->set_value(0,-double(Sign));
// Finish problem setup: Set the sign for the source function
// in all elements
//Find number of elements in mesh
unsigned n_element = mesh_pt()->nelement();
// Loop over the elements
for(unsigned i=0;i(mesh_pt()->element_pt(i));
// Set the sign of the source function
specific_element_pt->sign() = Sign;
}
//Assign the global and local equations numbers for the problem
cout << "Number of equations is " << assign_eqn_numbers() << std::endl;
}
/// Check that everything has been set up properly
void actions_before_newton_solve()
{
if (0==self_test())
{
cout << "Problem has been set up correctly and can be solved."
<< std::endl;
}
else
{
throw
OomphLibError("Trouble! Check error messages and fix the problems.\n",
"DemoPoissonProblem::actions_before_newton_solve()",
OOMPH_EXCEPTION_LOCATION);
}
}
/// Print out the result after the solve
void actions_after_newton_solve()
{
ofstream file ("result.dat");
mesh_pt()->output(file);
}
private:
/// The sign of the source function
int Sign;
}; //End of problem definition
//----------------------------MAIN FUNCTION-------------------------------
int main()
{
//Build the problem
DemoPoissonProblem problem;
//Solve the problem, using Newton's method
problem.newton_solve();
}