The boundary value problem for the first biharmonic equation is defined as follows:
\begin{eqnarray*}
\Delta^2 u &=& f \quad \text{in} \quad \Omega, \\
u &=& g_1 \quad \text{on} \quad \partial \Omega, \\
\mathbf{n} \cdot \nabla u &=& g_2 \quad \text{on} \quad \partial \Omega,
\end{eqnarray*}
where \(\mathbf{n}\) is the outward pointing unit normal vector. The boundary value problem for the second biharmonic equation is defined as follows:
\begin{eqnarray*}
\Delta^2 u &=& f \quad \text{in} \quad \Omega, \\
u &=& g_1 \quad \text{on} \quad \partial \Omega, \\
\Delta u &=& g_2 \quad \text{on} \quad \partial \Omega.
\end{eqnarray*}
The system to assemble a system of equations based on this weak formulation is constructed in biharmonic2_example.cpp using the gsExprAssembler. The biharmonic problem is only solvable on domains with a single patch.
Then we define our command line options. For example, we use the option -f
or --file
to set the path to the file that contains our geometry. If we wish to use the second biharmonic problem, then we use the option --second
to run the biharmonic problem with \(u=g_1\) and \(\Delta u=g_2\)
bool plot = false;
bool last = false;
bool second = false;
std::string fn;
gsCmdLine cmd("Example for solving the biharmonic problem (single patch only).");
cmd.addInt("p", "degree","Set discrete polynomial degree", degree);
cmd.addInt("s", "smoothness", "Set discrete regularity", smoothness);
cmd.addInt("l", "refinementLoop", "Number of refinement steps", numRefine);
cmd.addString("f", "file", "Input geometry file (with .xml)", fn);
cmd.addSwitch("last", "Solve problem only on the last level of h-refinement", last);
cmd.addSwitch("plot", "Create a ParaView visualization file with the solution", plot);
cmd.addSwitch("second", "Compute second biharmonic problem with u = g1 and Delta u = g2 "
"(default first biharmonic problem: u = g1 and partial_n u = g2)", second);
try { cmd.getValues(argc,argv); } catch (int rv) { return rv; }
#define index_t
Definition gsConfig.h:32
Then we read the geometry: Either it can be from a file that contains the geometry or we use the default one.
gsMultiPatch<> mp;
if (fn.empty())
else
{
gsInfo <<
"Filedata: " << fn <<
"\n";
gsReadFile<>(fn, mp);
}
mp.clearTopology();
mp.computeTopology();
if (mp.nPatches() != 1)
{
gsInfo <<
"The geometry has more than one patch. Run the code with a single patch!\n";
return EXIT_FAILURE;
}
#define gsInfo
Definition gsDebug.h:43
static TensorBSpline2Ptr BSplineFatQuarterAnnulus(T const &r0=1, T const &r1=2)
Definition gsNurbsCreator.hpp:949
Before we run the problem, we set the degree and/or refine the basis. This is done in the following snippet:
gsMultiBasis<> basis(mp, false);
basis.setDegree(degree);
if (last)
{
for (
index_t r =0; r < numRefine; ++r)
basis.uniformRefine(1, degree-smoothness);
numRefine = 0;
}
Then we set the boundary conditions. Depending on which biharmonic problem we choose, each boundary will be is set to the corresponding functions.
gsBoundaryConditions<> bc;
for (gsMultiPatch<>::const_biterator bit = mp.bBegin(); bit != mp.bEnd(); ++bit)
{
gsFunctionExpr<> laplace ("-16*pi*pi*(2*cos(4*pi*x)*cos(4*pi*y) - cos(4*pi*x) - cos(4*pi*y))",2);
gsFunctionExpr<> sol1der("-4*pi*(cos(4*pi*y) - 1)*sin(4*pi*x)",
"-4*pi*(cos(4*pi*x) - 1)*sin(4*pi*y)", 2);
if (second)
else
}
bc.setGeoMap(mp);
gsInfo <<
"Boundary conditions:\n" << bc <<
"\n";
@ dirichlet
Dirichlet type.
Definition gsBoundaryConditions.h:31
@ neumann
Neumann type.
Definition gsBoundaryConditions.h:33
@ laplace
Laplace type, e.g. \Delta u = g.
Definition gsBoundaryConditions.h:38
Now it is time to look at the setup of our boundary value problem. In the code snippet below, we first define our gsExprAssembler. The inputs (1,1) here mean that we have 1 test function space and 1 trial function space. Furthermore, few typedef
s can be seen. These are just there to prevent lengthy types. The assembler is fed with the basis using the gsExprAssembler::setIntegrationElements() function. When the basis is defined, a gsExprEvaluator is constructed for later use. The geometry is connected to the gsExprAssembler via the gsExprAssembler::getMap() function and we create an expression for the space via gsExprAssembler::getSpace(). Functions for the right-hand side and the manufactured solution are constructed via gsExprAssembler::getCoef() and gsExprEvaluator::getVariable() functions, respectively. Note that the manufactured solution is defined in the gsExprEvaluator because it will not be used in assembly of the system. The problem setup finishes by defining the solution object, which is constructed via gsExprAssembler::getSolution() which requires a space and a reference to a vector that will contain the solution coefficients later.
gsExprAssembler<real_t> A(1,1);
A.setIntegrationElements(basis);
gsExprEvaluator<real_t> ev(A);
auto G = A.getMap(mp);
auto ff = A.getCoeff(f, G);
auto u = A.getSpace(basis);
gsMatrix<real_t> solVector;
auto u_sol = A.getSolution(u, solVector);
auto u_ex = ev.getVariable(ms, G);
The next step is to actually solve the system at hand. In this particular example, this is done within a loop of uniform refinement (unless the last
option is selected) in order to study convergence behaviour of the method.
As seen in the snippet below, we first define the solver that we will use (based on the Eigen-library) and we define the vectors to store the errors.
Then, within the loop, we first refine the basis uniformly.
For the biharmonic problem, we fix the dofs in the function setMapperForBiharmonic(), which is defined in the beginning of the file (see the full file at the end).
Then the gsExprAssembler<>::space::setupMapper() function is used. This function initialize the matrix and the rhs depending on the mapper.
Then we compute the eliminated boundary values with the L2-projection where we call setDirichletNeumannValuesL2Projection(). Again this function is defined at the beginning of the file.
Then the gsExprAssembler is initialized using gsExprAssembler::initSystem() and using gsExprAssembler::assemble() the matrix expression (first argument) and the rhs expression (second argument) are assembled to a system.
When the System is assembled, the system is factorized using gsSparseSolver::compute() and solved via gsSparseSolver::solve(). The errors are then computed via the gsExprEvaluator, where it should be noted that the solution object is appearing, which uses the solution vector and the space as defined earlier.
gsSparseSolver<real_t>::SimplicialLDLT solver;
gsVector<real_t> l2err(numRefine+1), h1err(numRefine+1), h2err(numRefine+1),
dofs(numRefine+1), meshsize(numRefine+1);
gsInfo<<
"(dot1=assembled, dot2=solved, dot3=got_error)\n"
"\nDoFs: ";
double setup_time(0), ma_time(0), slv_time(0), err_time(0);
for (
index_t r=0; r<=numRefine; ++r)
{
basis.uniformRefine(1,degree -smoothness);
meshsize[r] = basis.basis(0).getMaxCellLength();
gsDofMapper map;
setMapperForBiharmonic(bc, basis,map);
u.setupMapper(map);
setDirichletNeumannValuesL2Projection(mp, basis, bc, u);
A.initSystem();
setup_time += timer.stop();
gsInfo<< A.numDofs() <<std::flush;
timer.restart();
A.assemble(ilapl(u, G) * ilapl(u, G).tr() * meas(G),u * ff * meas(G));
auto g_L = A.getBdrFunction(G);
A.assembleBdr(bc.get("Laplace"), (igrad(u, G) * nv(G)) * g_L.tr() );
dofs[r] = A.numDofs();
ma_time += timer.stop();
timer.restart();
solver.compute( A.matrix() );
solVector = solver.solve(A.rhs());
slv_time += timer.stop();
timer.restart();
l2err[r]= math::sqrt( ev.integral( (u_ex - u_sol).sqNorm() * meas(G) ) );
h1err[r]= l2err[r] +
math::sqrt(ev.integral( ( igrad(u_ex) - igrad(u_sol,G) ).sqNorm() * meas(G) ));
h2err[r]= h1err[r] +
math::sqrt(ev.integral( ( ihess(u_ex) - ihess(u_sol,G) ).sqNorm() * meas(G) ));
err_time += timer.stop();
}
gsGenericStopwatch< WallClock > gsStopwatch
A stop-watch measuring real (wall) time.
Definition gsStopwatch.h:160
When the loop is runnig several times, we can compute the L2-error, H1-error and H2-error. The expected rates should be p+1 for the L2-error, p for the H1-error and p-1 for the H2-error.
gsInfo <<
"\nL2 error: "<<std::scientific<<std::setprecision(3)<<l2err.transpose()<<
"\n";
gsInfo <<
"H1 error: "<<std::scientific<<h1err.transpose()<<
"\n";
gsInfo <<
"H2 error: "<<std::scientific<<h2err.transpose()<<
"\n";
if (numRefine>0)
{
gsInfo<<
"\nEoC (L2): " << std::fixed<<std::setprecision(2)
<< ( l2err.head(numRefine).array() /
l2err.tail(numRefine).array() ).log().transpose() / std::log(2.0)
<<"\n";
gsInfo<<
"EoC (H1): "<< std::fixed<<std::setprecision(2)
<<( h1err.head(numRefine).array() /
h1err.tail(numRefine).array() ).log().transpose() / std::log(2.0) <<"\n";
gsInfo<<
"EoC (H2): "<< std::fixed<<std::setprecision(2)
<<( h2err.head(numRefine).array() /
h2err.tail(numRefine).array() ).log().transpose() / std::log(2.0) <<"\n";
}
If we set the flag --plot
then the export to ParaView is done in the snipped below. The output file is solution.pvd.
if (plot)
{
gsInfo <<
"Plotting in Paraview...\n";
ev.options().setSwitch("plot.elements", true);
ev.writeParaview( u_sol , G, "solution");
gsInfo <<
"Saved with solution.pvd \n";
}
else
gsInfo <<
"Done. No output created, re-run with --plot to get a ParaView "
"file containing the solution.\n";
Annotated source file
Here is the full file examples/biharmonic2_example.cpp
. Clicking on a function or class name will lead you to its reference documentation.
{
for (gsBoxTopology::const_iiterator it = basis.topology().iBegin();
it != basis.topology().iEnd(); ++it)
{
basis.matchInterface(*it, mapper);
}
for (typename gsBoundaryConditions<real_t>::const_iterator
it = bc.
begin(
"Dirichlet"); it != bc.
end(
"Dirichlet"); ++it)
{
bnd = basis.basis(it->ps.patch).boundary(it->ps.side());
}
for (typename gsBoundaryConditions<real_t>::const_iterator
it = bc.
begin(
"Neumann"); it != bc.
end(
"Neumann"); ++it)
{
bnd = basis.basis(it->ps.patch).boundaryOffset(it->ps.side(),1);
}
}
{
for (gsBoxTopology::const_iiterator it = basis.topology().iBegin();
it != basis.topology().iEnd(); ++it)
{
basis.matchInterface(*it, mapperBdy);
}
for (size_t np = 0; np < mp.nPatches(); np++)
{
mapperBdy.markBoundary(np, bnd, 0);
}
mapperBdy.finalize();
A.setIntegrationElements(basis);
auto G = A.getMap(mp);
auto uu = A.getSpace(basis);
auto g_bdy = A.getBdrFunction(G);
uu.setupMapper(mapperBdy);
fixedDofs_A.setZero( uu.mapper().boundarySize(), 1 );
real_t lambda = 1e-5;
A.initSystem();
A.assembleBdr(bc.
get(
"Dirichlet"), uu * uu.tr() *
meas(G));
A.assembleBdr(bc.
get(
"Dirichlet"), uu * g_bdy *
meas(G));
A.assembleBdr(bc.
get(
"Neumann"),
lambda * (igrad(uu, G) *
nv(G).normalized()) * (igrad(uu, G) *
nv(G).normalized()).tr() *
meas(G));
A.assembleBdr(bc.
get(
"Neumann"),
lambda * (igrad(uu, G) *
nv(G).normalized()) * (g_bdy.tr() *
nv(G).normalized()) *
meas(G));
gsSparseSolver<real_t>::SimplicialLDLT solver;
solver.compute( A.matrix() );
for (size_t np = 0; np < mp.nPatches(); np++)
{
bnd.array() += sz;
for (
index_t i = 0; i < bnd.rows(); i++)
{
index_t ii = mapperBdy.asVector()(bnd(i,0));
}
sz += mapperBdy.patchSize(np,0);
}
}
int main(int argc, char *argv[])
{
bool plot = false;
bool last = false;
bool second = false;
std::string fn;
gsCmdLine cmd(
"Example for solving the biharmonic problem (single patch only).");
cmd.addInt("p", "degree","Set discrete polynomial degree", degree);
cmd.addInt("s", "smoothness", "Set discrete regularity", smoothness);
cmd.addInt("l", "refinementLoop", "Number of refinement steps", numRefine);
cmd.addString("f", "file", "Input geometry file (with .xml)", fn);
cmd.addSwitch("last", "Solve problem only on the last level of h-refinement", last);
cmd.addSwitch("plot", "Create a ParaView visualization file with the solution", plot);
cmd.addSwitch("second", "Compute second biharmonic problem with u = g1 and Delta u = g2 "
"(default first biharmonic problem: u = g1 and partial_n u = g2)", second);
try { cmd.getValues(argc,argv); } catch (int rv) { return rv; }
if (fn.empty())
else
{
gsInfo <<
"Filedata: " << fn <<
"\n";
}
mp.clearTopology();
mp.computeTopology();
if (mp.nPatches() != 1)
{
gsInfo <<
"The geometry has more than one patch. Run the code with a single patch!\n";
return EXIT_FAILURE;
}
gsFunctionExpr<>f(
"256*pi*pi*pi*pi*(4*cos(4*pi*x)*cos(4*pi*y) - cos(4*pi*x) - cos(4*pi*y))",2);
gsInfo <<
"Source function: " << f <<
"\n";
gsInfo <<
"Exact function: " << ms <<
"\n";
basis.setDegree(degree);
if (last)
{
for (
index_t r =0; r < numRefine; ++r)
basis.uniformRefine(1, degree-smoothness);
numRefine = 0;
}
for (gsMultiPatch<>::const_biterator bit = mp.bBegin(); bit != mp.bEnd(); ++bit)
{
gsFunctionExpr<> laplace (
"-16*pi*pi*(2*cos(4*pi*x)*cos(4*pi*y) - cos(4*pi*x) - cos(4*pi*y))",2);
"-4*pi*(cos(4*pi*x) - 1)*sin(4*pi*y)", 2);
if (second)
else
}
gsInfo <<
"Boundary conditions:\n" << bc <<
"\n";
A.setIntegrationElements(basis);
auto G = A.getMap(mp);
auto ff = A.getCoeff(f, G);
auto u = A.getSpace(basis);
auto u_sol = A.getSolution(u, solVector);
auto u_ex = ev.getVariable(ms, G);
#ifdef _OPENMP
gsInfo <<
"Available threads: "<< omp_get_max_threads() <<
"\n";
#endif
gsSparseSolver<real_t>::SimplicialLDLT solver;
dofs(numRefine+1), meshsize(numRefine+1);
gsInfo<<
"(dot1=assembled, dot2=solved, dot3=got_error)\n"
"\nDoFs: ";
double setup_time(0), ma_time(0), slv_time(0), err_time(0);
for (
index_t r=0; r<=numRefine; ++r)
{
basis.uniformRefine(1,degree -smoothness);
meshsize[r] = basis.basis(0).getMaxCellLength();
setMapperForBiharmonic(bc, basis,map);
u.setupMapper(map);
setDirichletNeumannValuesL2Projection(mp, basis, bc, u);
A.initSystem();
setup_time += timer.
stop();
gsInfo<< A.numDofs() <<std::flush;
A.assemble(ilapl(u, G) * ilapl(u, G).tr() *
meas(G),u * ff *
meas(G));
auto g_L = A.getBdrFunction(G);
A.assembleBdr(bc.get(
"Laplace"), (igrad(u, G) *
nv(G)) * g_L.tr() );
dofs[r] = A.numDofs();
solver.compute( A.matrix() );
solVector = solver.solve(A.rhs());
slv_time += timer.
stop();
l2err[r]= math::sqrt( ev.integral( (u_ex - u_sol).sqNorm() *
meas(G) ) );
h1err[r]= l2err[r] +
math::sqrt(ev.integral( ( igrad(u_ex) - igrad(u_sol,G) ).sqNorm() *
meas(G) ));
h2err[r]= h1err[r] +
math::sqrt(ev.integral( ( ihess(u_ex) - ihess(u_sol,G) ).sqNorm() *
meas(G) ));
err_time += timer.
stop();
}
gsInfo <<
"\n\nTotal time: "<< setup_time+ma_time+slv_time+err_time <<
"\n";
gsInfo <<
" Setup: "<< setup_time <<
"\n";
gsInfo <<
" Assembly: "<< ma_time <<
"\n";
gsInfo <<
" Solving: "<< slv_time <<
"\n";
gsInfo <<
" Norms: "<< err_time <<
"\n";
gsInfo <<
" Mesh-size: "<< meshsize.transpose() <<
"\n";
gsInfo <<
" Dofs: "<<dofs.transpose() <<
"\n";
gsInfo <<
"\nL2 error: "<<std::scientific<<std::setprecision(3)<<l2err.transpose()<<
"\n";
gsInfo <<
"H1 error: "<<std::scientific<<h1err.transpose()<<
"\n";
gsInfo <<
"H2 error: "<<std::scientific<<h2err.transpose()<<
"\n";
if (numRefine>0)
{
gsInfo<<
"\nEoC (L2): " << std::fixed<<std::setprecision(2)
<< ( l2err.head(numRefine).array() /
l2err.tail(numRefine).array() ).log().transpose() / std::log(2.0)
<<"\n";
gsInfo<<
"EoC (H1): "<< std::fixed<<std::setprecision(2)
<<( h1err.head(numRefine).array() /
h1err.tail(numRefine).array() ).log().transpose() / std::log(2.0) <<"\n";
gsInfo<<
"EoC (H2): "<< std::fixed<<std::setprecision(2)
<<( h2err.head(numRefine).array() /
h2err.tail(numRefine).array() ).log().transpose() / std::log(2.0) <<"\n";
}
if (plot)
{
gsInfo <<
"Plotting in Paraview...\n";
ev.options().setSwitch("plot.elements", true);
ev.writeParaview( u_sol , G, "solution");
gsInfo <<
"Saved with solution.pvd \n";
}
else
gsInfo <<
"Done. No output created, re-run with --plot to get a ParaView "
"file containing the solution.\n";
return EXIT_SUCCESS;
}
Definition gsExpressions.h:973
Class containing a set of boundary conditions.
Definition gsBoundaryConditions.h:342
void addCondition(int p, boxSide s, condition_type::type t, gsFunctionSet< T > *f, short_t unknown=0, bool parametric=false, int comp=-1)
Adds another boundary condition.
Definition gsBoundaryConditions.h:650
bcRefList get(const std::string &label, const short_t unk=0, int comp=-1) const
Definition gsBoundaryConditions.h:420
const_iterator begin(const std::string &label) const
Returns a const-iterator to the beginning of the Bc container of type label.
Definition gsBoundaryConditions.h:471
const_iterator end(const std::string &label) const
Returns a const-iterator to the end of the Bc container of type label.
Definition gsBoundaryConditions.h:477
void setGeoMap(const gsFunctionSet< T > &gm)
Set the geometry map to evaluate boundary conditions.
Definition gsBoundaryConditions.h:916
Class for command-line argument parsing.
Definition gsCmdLine.h:57
Maintains a mapping from patch-local dofs to global dof indices and allows the elimination of individ...
Definition gsDofMapper.h:69
void markBoundary(index_t k, const gsMatrix< index_t > &boundaryDofs, index_t comp=0)
Mark the local dofs boundaryDofs of patch k as eliminated.
Definition gsDofMapper.cpp:188
gsVector< index_t > findFree(const index_t k) const
Returns all free dofs on patch k (local dof indices)
Definition gsDofMapper.cpp:652
void init(const gsMultiBasis< T > &bases, index_t nComp=1)
Initialize by a gsMultiBasis.
Definition gsDofMapper.hpp:20
void finalize()
Must be called after all boundaries and interfaces have been marked to set up the dof numbering.
Definition gsDofMapper.cpp:240
index_t boundarySize() const
Returns the number of eliminated dofs.
Definition gsDofMapper.h:457
index_t global_to_bindex(index_t gl) const
Returns the boundary index of global dof gl.
Definition gsDofMapper.h:364
gsVector< index_t > asVector(index_t comp=0) const
Returns a vector taking flat local indices to global.
Definition gsDofMapper.cpp:80
Definition gsExprAssembler.h:33
Generic evaluator of isogeometric expressions.
Definition gsExprEvaluator.h:39
Class defining a multivariate (real or vector) function given by a string mathematical expression.
Definition gsFunctionExpr.h:52
A Stopwatch object can be used to measure execution time of code, algorithms, etc.
Definition gsStopwatch.h:73
double stop()
Return elapsed time in seconds.
Definition gsStopwatch.h:83
void restart()
Start taking the time.
Definition gsStopwatch.h:80
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
Holds a set of patch-wise bases and their topology information.
Definition gsMultiBasis.h:37
Container class for a set of geometry patches and their topology, that is, the interface connections ...
Definition gsMultiPatch.h:100
Reads an object from a data file, if such the requested object exists in the file.
Definition gsReadFile.h:43
A vector with arbitrary coefficient type and fixed or dynamic size.
Definition gsVector.h:37
Main header to be included by clients using the G+Smo library.
EIGEN_STRONG_INLINE meas_expr< T > meas(const gsGeometryMap< T > &G)
The measure of a geometry map.
Definition gsExpressions.h:4555
EIGEN_STRONG_INLINE onormal_expr< T > nv(const gsGeometryMap< T > &u)
The (outer pointing) boundary normal of a geometry map.
Definition gsExpressions.h:4506
The G+Smo namespace, containing all definitions for the library.
Class gsNurbsCreator provides some simple examples of Nurbs Geometries.
Definition gsNurbsCreator.h:37