Hello Deal II community,
I have just started with deal ii, so I am trying to solve a basic laplace
equation on a square domain with dirichilet bc on two faces and neumann bc
on the remaining two. I have used steps from step-7 to implement neumann
bc, but i am running into the following error:
--------------------------------------------------------
An error occurred in line <133> of file
</home/umair/dealdotii/dealii-9.1.1/include/deal.II/base/function.templates.h>
in function
dealii::Tensor<1, dim, Number> dealii::Function<dim,
RangeNumberType>::gradient(const dealii::Point<dim>&, unsigned int) const
[with int dim = 2; RangeNumberType = double]
The violated condition was:
false
Additional information:
You (or a place in the library) are trying to call a function that is
declared as a virtual function in a base class but that has not been
overridden in your derived class.
This exception happens in cases where the base class cannot provide a
useful default implementation for the virtual function, but where we also
do not want to mark the function as abstract (i.e., with '=0' at the end)
because the function is not essential to the class in many contexts. In
cases like this, the base class provides a dummy implementation that makes
the compiler happy, but that then throws the current exception.
A concrete example would be the 'Function' class. It declares the existence
of 'value()' and 'gradient()' member functions, and both are marked as
'virtual'. Derived classes have to override these functions for the values
and gradients of a particular function. On the other hand, not every
function has a gradient, and even for those that do, not every program
actually needs to evaluate it. Consequently, there is no *requirement* that
a derived class actually override the 'gradient()' function (as there would
be had it been marked as abstract). But, since the base class cannot know
how to compute the gradient, if a derived class does not override the
'gradient()' function and it is called anyway, then the default
implementation in the base class will simply throw an exception.
The exception you see is what happens in cases such as the one just
illustrated. To fix the problem, you need to investigate whether the
function being called should indeed have been called; if the answer is
'yes', then you need to implement the missing override in your class.
Stacktrace:
-----------
#0 /usr/local/lib/libdeal_II.g.so.9.1.1: dealii::Function<2,
double>::gradient(dealii::Point<2, double> const&, unsigned int) const
#1 /home/umair/CLionProjects/MyLaplace2/cmake-build-debug/MyLaplace2:
laplace_solver::assemble_system()
#2 /home/umair/CLionProjects/MyLaplace2/cmake-build-debug/MyLaplace2:
laplace_solver::run()
#3 /home/umair/CLionProjects/MyLaplace2/cmake-build-debug/MyLaplace2: main
--------------------------------------------------------
I have attached my code for more help.
Thank you in advance,
Umair.
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/2fc8eb10-f1e3-4706-80a1-9837450f7874n%40googlegroups.com.
//
// Created by umair on 08/07/20.
//
#ifndef UNTITLED_LAPLACE_SOLVER_H
#define UNTITLED_LAPLACE_SOLVER_H
#include <deal.II/grid/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/numerics/data_out.h>
#include <fstream>
#include <iostream>
#include "RightHandSide.h"
#include "BoundaryValues.h"
#include "Neumann.h"
using namespace dealii;
class laplace_solver {
public:
laplace_solver ();
void run ();
private:
void make_grid ();
void setup_system ();
void assemble_system ();
void solve ();
void output_results () const;
Triangulation<2> triangulation;
FE_Q<2> fe;
DoFHandler<2> dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};
laplace_solver::laplace_solver()
: fe(1) // linear elements
, dof_handler(triangulation)
{}
void laplace_solver::make_grid()
{
GridGenerator::hyper_cube(triangulation, 0, 1);
triangulation.refine_global(2);
std::cout << "Number of active cells: " << triangulation.n_active_cells()
<< std::endl;
for (auto &cell : triangulation.active_cell_iterators())
for (unsigned int face_n = 0;
face_n < GeometryInfo<2>::faces_per_cell;
++face_n)
if (cell->at_boundary(face_n)){
if (std::fabs(cell->face(face_n)->center()[0] - 1) < 1e-10 || std::fabs(cell->face(face_n)->center()[1] - 1) < 1e-10)
cell->face(face_n)->set_boundary_id(55);
else if (std::fabs(cell->face(face_n)->center()[0] - 0) < 1e-10 || std::fabs(cell->face(face_n)->center()[1] - 0) < 1e-10)
cell->face(face_n)->set_boundary_id(555);
}
triangulation.refine_global(1);
}
void laplace_solver::setup_system()
{
dof_handler.distribute_dofs(fe);
std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
DynamicSparsityPattern dsp(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler, dsp);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
}
void laplace_solver::assemble_system()
{
QGauss<2> quadrature_formula(fe.degree + 1);
QGauss<1> face_quadrature_formula(fe.degree + 1); //For Neumann BC
RightHandSide right_hand_side;
FEValues<2> fe_values(fe,
quadrature_formula,
update_values | update_gradients | update_quadrature_points | update_JxW_values);
FEFaceValues<2> fe_face_values(fe,
face_quadrature_formula,
update_values | update_quadrature_points |
update_normal_vectors |
update_JxW_values); //For Neumann BC
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
const unsigned int n_face_q_points = face_quadrature_formula.size(); //For Neumann BC
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
Neumann neumann_bc; //initialising Function for Neumann BC
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
cell_matrix = 0;
cell_rhs = 0;
for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < dofs_per_cell; ++j)
cell_matrix(i, j) +=
(fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
fe_values.JxW(q_index)); // dx
const auto x_q = fe_values.quadrature_point(q_index);
cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
right_hand_side.value(x_q) * // f(x_q)
fe_values.JxW(q_index)); // dx
}
// Adding Boundary integral term for Neumann BC
for (unsigned int face_number = 0;
face_number < GeometryInfo<2>::faces_per_cell;
++face_number)
if (cell->face(face_number)->at_boundary() &&
(cell->face(face_number)->boundary_id() == 555))
{
fe_face_values.reinit(cell, face_number);
for (unsigned int q_point = 0; q_point < n_face_q_points;
++q_point)
{
const double neumann_value =
(neumann_bc.gradient(
fe_face_values.quadrature_point(q_point)) *
fe_face_values.normal_vector(q_point));
for (unsigned int i = 0; i < dofs_per_cell; ++i)
cell_rhs(i) +=
(neumann_value * // g(x_q)
fe_face_values.shape_value(i, q_point) * // phi_i(x_q)
fe_face_values.JxW(q_point)); // dx
}
}
cell->get_dof_indices(local_dof_indices);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
system_matrix.add(local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i, j));
for (unsigned int i = 0; i < dofs_per_cell; ++i)
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
std::map<types::global_dof_index, double> boundary_values;
VectorTools::interpolate_boundary_values(dof_handler,
55,
BoundaryValues(),
boundary_values);
MatrixTools::apply_boundary_values(boundary_values,
system_matrix,
solution,
system_rhs);
}
void laplace_solver::solve()
{
SolverControl solver_control(1000, 1e-12);
SolverCG<Vector<double>> solver(solver_control);
solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
}
void laplace_solver::output_results() const
{
DataOut<2> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "solution");
data_out.build_patches();
std::ofstream output("solution.vtk");
data_out.write_vtk(output);
}
void laplace_solver::run()
{
std::cout<<"The RUN METHOD is called"<< std::endl;
make_grid();
setup_system();
assemble_system();
solve();
output_results();
}
#endif //UNTITLED_LAPLACE_SOLVER_H
//
// Created by umair on 21/08/20.
//
#ifndef MYLAPLACE2_NEUMANN_H
#define MYLAPLACE2_NEUMANN_H
#include <deal.II/base/function.h>
#include <cmath>
using namespace dealii;
class Neumann : public Function<2>
{
public:
virtual double value(const Point<2> & p,
const unsigned int component = 0) const override;
};
double Neumann::value(const Point<2> &p,
const unsigned int /*component*/) const {
return std::sin(5*p(0));
}
#endif //MYLAPLACE2_NEUMANN_H
cmake_minimum_required(VERSION 3.16)
# Set the name of the project and target:
SET(TARGET "MyLaplace2")
#set(CMAKE_CXX_STANDARD 14)
SET(TARGET_SRC
main.cpp
)
FIND_PACKAGE(deal.II 9.1.0 QUIET
HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
)
IF(NOT ${deal.II_FOUND})
MESSAGE(FATAL_ERROR "\n"
"*** Could not locate a (sufficiently recent) version of deal.II.
***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II
to cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this
path."
)
ENDIF()
DEAL_II_INITIALIZE_CACHED_VARIABLES()
PROJECT(${TARGET})
DEAL_II_INVOKE_AUTOPILOT()#include <iostream>
#include "laplace_solver.h"
int main() {
std::cout << "Hello, World!" << std::endl;
laplace_solver laplace_problem;
laplace_problem.run();
return 0;
// Comment need to be stripped
}
//
// Created by umair on 13/07/20.
//
#ifndef MYLAPLACE_BOUNDARYVALUES_H
#define MYLAPLACE_BOUNDARYVALUES_H
#include <deal.II/base/function.h>
using namespace dealii;
class BoundaryValues : public Function<2>
{
public:
virtual double value(const Point<2> & p,
const unsigned int component = 0) const override;
};
double BoundaryValues::value(const Point<2> &p,
const unsigned int /*component*/) const
{
return (0.0);
}
#endif //MYLAPLACE_BOUNDARYVALUES_H
//
// Created by umair on 13/07/20.
//
#ifndef MYLAPLACE_RIGHTHANDSIDE_H
#define MYLAPLACE_RIGHTHANDSIDE_H
#include <deal.II/base/function.h>
using namespace dealii;
class RightHandSide : public Function<2>
{
public:
virtual double value(const Point<2> & p,
const unsigned int component = 0) const override;
};
double RightHandSide::value(const Point<2> &p,
const unsigned int /*component*/) const {
return 1.0;
}
#endif //MYLAPLACE_RIGHTHANDSIDE_H
