Benjamin Kirk wrote:
> I don't seem to see the attachment... But I have to ask:
>
> If the assembly function was taken from example 3 (2 does no assembly(?)),
> are you constraining the hanging degrees of freedom? That is, is there a
> line like
>
> // We have now built the element matrix and RHS vector in terms
> // of the element degrees of freedom. However, it is possible
> // that some of the element DOFs are constrained to enforce
> // solution continuity, i.e. they are not really "free". We need
> // to constrain those DOFs in terms of non-constrained DOFs to
> // ensure a continuous solution. The
> // \p DofMap::constrain_element_matrix_and_vector() method does
> // just that.
> dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
>
> (c.f. example 10)
>
> -Ben
>
It is example 3; since the attachment didn't pass through i give the
whole program below.
I tried also constraining element matrices but it doesn't make any
difference.
Here is the program:
Mladen
#include <iostream>
#include <algorithm>
#include <cmath>
#include <cassert>
#include <cstdlib>
#include "libmesh.h"
#include "mesh.h"
#include "mesh_generation.h"
#include "mesh_refinement.h"
#include "linear_implicit_system.h"
#include "equation_systems.h"
#include "vector_value.h"
#include "exact_solution.h"
#include "fe.h"
#include "quadrature_gauss.h"
#include "dof_map.h"
#include "sparse_matrix.h"
#include "numeric_vector.h"
#include "dense_matrix.h"
#include "dense_vector.h"
#include "elem.h"
#include "string_to_enum.h"
void assemble_poisson(EquationSystems& es, const std::string& system_name);
Number es_ptr(const Point &p, const Parameters &Parameters,
const std::string &sys_name, const std::string &unknown_name)
{
const Real pi = libMesh::pi;
double x=p(0), y=p(1);
return std::sin( pi * x) * std::sin( pi * y);
}
Gradient ges_ptr(const Point &p, const Parameters ¶meters,
const std::string &sys_name, const std::string &unknown_name)
{
const Real pi = libMesh::pi;
double x=p(0), y=p(1);
Gradient g;
g(0) = pi * std::cos( pi * x) * std::sin( pi * y);
g(1) = pi * std::sin( pi * x) * std::cos( pi * y);
g(2) = 0.0;
return g;
}
//==============================================================
Real exact(double x, double y, double z)
{
const Real pi = libMesh::pi;
return std::sin( pi * x) * std::sin( pi * y);
}
double exact(Point p)
{
const Real pi = libMesh::pi;
double x=p(0), y=p(1);
return std::sin( pi * x) * std::sin( pi * y);
}
//==============================================================
int main (int argc, char** argv)
{
using std::cout;
using std::endl;
LibMeshInit init (argc, argv);
Mesh mesh (2);
ElemType el_type = Utility::string_to_enum<ElemType>("TRI3");
unsigned int n_refinement = 4;
MeshTools::Generation::build_square (mesh,
10, 10,
0.0, 1.0,
0.0, 1.0,
el_type);
MeshRefinement mesh_refinement(mesh);
EquationSystems es(mesh);
LinearImplicitSystem & system =
es.add_system<LinearImplicitSystem>("Poisson");
system.add_variable ("p", Utility::string_to_enum<Order> ("FIRST"));
system.attach_assemble_function (assemble_poisson);
es.parameters.set<Real>("linear solver tolerance") =1.0E-15;
es.init();
ExactSolution ex_sol(es);
ex_sol.attach_exact_value(es_ptr);
ex_sol.attach_exact_deriv(ges_ptr);
ex_sol.extra_quadrature_order(5);
Real l2_error=0.0, l2_error_old=0.0;
Real linf_error=0.0, linf_error_old=0.0;
Real h1_error=0.0, h1_error_old=0.0;
std::ofstream out_file("results.txt"); // output file
// out_file.open("results.txt");
out_file << " E R R O R C O N V E R G E
N C E R A T E \n";
out_file << " H1 L2 Linf
H1 L2 Linf \n";
out_file << "========== ========= ========== =========
========= ==========\n";
int max_refinement=n_refinement;
for(int refinement=0; refinement < max_refinement; ++ refinement)
{
system.solve();
ex_sol.compute_error("Poisson","p");
l2_error = ex_sol.l2_error("Poisson","p");
h1_error = ex_sol.h1_error("Poisson","p");
linf_error = ex_sol.l_inf_error("Poisson","p");
out_file << std::scientific
<< h1_error << " " << l2_error << " "<< linf_error<< " ";
if(refinement >0)
out_file << std::fixed
<< std::log2(h1_error_old/h1_error)<< " "
<< std::log2(l2_error_old/l2_error) << " "
<< std::log2(linf_error_old/linf_error);
out_file <<std::endl;
h1_error_old = h1_error;
l2_error_old = l2_error;
linf_error_old = linf_error;
mesh_refinement.uniformly_refine();
es.reinit();
}
out_file.close();
return 0;
}
void assemble_poisson(EquationSystems& es, const std::string& system_name)
{
assert (system_name == "Poisson");
const MeshBase& mesh = es.get_mesh();
const unsigned int dim = mesh.mesh_dimension();
LinearImplicitSystem& system =
es.get_system<LinearImplicitSystem>("Poisson");
const DofMap& dof_map = system.get_dof_map();
FEType fe_type = dof_map.variable_type(0); // 0 = var_number
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
QGauss qrule (dim, FIFTH);
fe->attach_quadrature_rule (&qrule);
AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
QGauss qface(dim-1, FIFTH);
fe_face->attach_quadrature_rule (&qface);
const std::vector<Real>& JxW = fe->get_JxW();
const std::vector<Point>& q_point = fe->get_xyz();
const std::vector<std::vector<Real> >& phi = fe->get_phi();
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
DenseMatrix<Number> Ke;
DenseVector<Number> Fe;
std::vector<unsigned int> dof_indices;
MeshBase::const_element_iterator el =
mesh.local_elements_begin();
const MeshBase::const_element_iterator end_el =
mesh.local_elements_end();
for ( ; el != end_el; ++el)
{
const Elem* elem = *el;
dof_map.dof_indices (elem, dof_indices);
unsigned int n_dofs = dof_indices.size();
fe->reinit (elem);
Ke.resize (n_dofs, n_dofs);
Fe.resize (n_dofs);
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]);
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
const Real x = q_point[qp](0);
const Real y = q_point[qp](1);
const Real z = q_point[qp](2);
const Real eps = 1.e-5;
const Real uxx = (exact(x-eps,y,z) + exact(x+eps,y,z)
-2.*exact(x,y,z))/eps/eps;
const Real uyy = (exact(x,y-eps,z) + exact(x,y+eps,z)
-2.*exact(x,y,z))/eps/eps;
Real fxy;
fxy = - (uxx + uyy);
for (unsigned int i=0; i<phi.size(); i++)
Fe(i) += JxW[qp]*fxy*phi[i][qp];
}
{
for (unsigned int side=0; side<elem->n_sides(); side++)
if (elem->neighbor(side) == NULL)
{
const Real penalty = 1.e20;
const std::vector<std::vector<Real> >& phi_face =
fe_face->get_phi();
const std::vector<Real>& JxW_face = fe_face->get_JxW();
const std::vector<Point >& qface_point = fe_face->get_xyz();
fe_face->reinit(elem, side);
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
const Real xf = qface_point[qp](0);
const Real yf = qface_point[qp](1);
const Real zf = qface_point[qp](2);
const Real value = exact(xf, yf, zf);
for (unsigned int i=0; i<phi_face.size(); i++)
for (unsigned int j=0; j<phi_face.size(); j++)
Ke(i,j) +=
JxW_face[qp]*penalty*phi_face[i][qp]*phi_face[j][qp];
for (unsigned int i=0; i<phi_face.size(); i++)
Fe(i) += JxW_face[qp]*penalty*value*phi_face[i][qp];
}
}
}
system.matrix->add_matrix (Ke, dof_indices);
system.rhs->add_vector (Fe, dof_indices);
}
}
-------------------------------------------------------------------------
This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
Build the coolest Linux based applications with Moblin SDK & win great prizes
Grand prize is a trip for two to an Open Source event anywhere in the world
http://moblin-contest.org/redirect.php?banner_id=100&url=/
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
