Hi everyone,
I am trying to implement a Dirichlet boundary condition for a problem very
similar to step-44. The boundary condition is just a displacement in the +x
direction on the +x face of a cube. However, when I call
VectorTools::interpolate_boundary_values, it seems it's calling a default
(zero?) version of Function<dim>::value rather than
IncrementalDisplacement<dim>::value, which should override the default
value function.
Any help on how to solve this issue would be appreciated.
Attached is a minimal working example. In my opinion, if my overridden
function is being called, then the cout statement inside of it ("Overridden
value function called") should be displayed, right?
Some things I've tried:
- Removing the "override" keyword and see if deal.II would raise an
exception. The program runs entirely and nothing is displayed on the
terminal.
- Calling the function IncrementalDisplacement<dim>::value directly.
This works fine.
- Overloading vector_value and vector_value_list. Same thing, these
functions are not being called either.
- Stepping into the VectorTools::interpolate_boundary_values function
using gdb, but the output is the following:
(gdb) step
dealii::VectorTools::interpolate_boundary_values<3, 3, double> (dof=...,
boundary_component=1, boundary_function=..., constraints=...,
component_mask=...)
at
/home/javier/dealii-9.3/tmp/unpack/deal.II-v9.3.2/include/deal.II/numerics/vector_tools_boundary.templates.h:603
603
/home/javier/dealii-9.3/tmp/unpack/deal.II-v9.3.2/include/deal.II/numerics/vector_tools_boundary.templates.h:
No such file or directory.
Thanks!
Javier
--
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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/b88e509a-8535-47f2-b49e-54ac8f1fa16bn%40googlegroups.com.
// Using deal.II 9.3.2
#include <deal.II/base/function.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_dgp_monomial.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <iostream>
namespace MainProblem
{
using namespace dealii;
template <int dim>
class IncrementalDisplacement : public Function<dim>
{
public:
IncrementalDisplacement(const double a, const double b)
:
Function<dim>(dim + 2),
a(a),
b(b)
{}
virtual double value (const Point<dim> &/*p*/,
const unsigned int component = 0) const override;
virtual void vector_value (const Point<dim> &/*p*/,
Vector<double> &values) const override;
virtual void vector_value_list (const std::vector<Point<dim>> &points,
std::vector<Vector<double>> &value_list) const override;
private:
const double a, b;
};
template <>
double IncrementalDisplacement<3>::value(const Point<3> &/*p*/,
const unsigned int component) const
{
std::cout << "Overridden value function called" << std::endl;
return (component == 0) ? a - b : 0.0;
}
template <>
void IncrementalDisplacement<3>::vector_value (const Point<3> &/*p*/,
Vector<double> &values) const
{
std::cout << "Overridden vector_value function called" << std::endl;
values = 0;
values(0) = a - b;
}
template <>
void IncrementalDisplacement<3>::vector_value_list (const std::vector<Point<3>> &points,
std::vector<Vector<double>> &value_list) const
{
std::cout << "Overridden vector_value_list function called" << std::endl;
const unsigned int n_points = points.size();
Assert (value_list.size() == n_points,
ExcDimensionMismatch (value_list.size(), n_points));
for (unsigned int p=0; p < n_points; ++p)
IncrementalDisplacement<3>::vector_value(points[p], value_list[p]);
}
}
int main()
{
using namespace MainProblem;
const int dim = 3;
// Generate unit cube for testing
Triangulation<dim> triangulation;
GridGenerator::hyper_cube(triangulation);
triangulation.refine_global(2);
// and a dummy dof_handler element
DoFHandler<dim> dof_handler(triangulation);
// The context is the problem solved in step-44
// (hyperelasticity using a three-field formulation)
const FESystem<dim> fe(FE_Q<dim>(2), dim, // displacement
FE_DGPMonomial<dim>(1), 1, // pressure
FE_DGPMonomial<dim>(1), 1); // dilation
dof_handler.distribute_dofs(fe);
// Add a dummy constraint object
AffineConstraints<double> constraints;
// This would be part of the make_constraints function in step-44
const FEValuesExtractors::Scalar x_displacement(0);
const FEValuesExtractors::Scalar y_displacement(1);
const FEValuesExtractors::Scalar z_displacement(2);
IncrementalDisplacement<dim> inc_object(3.0, 2.0);
// Want to apply boundary condition on the +x face of the cube
constraints.clear();
VectorTools::interpolate_boundary_values(dof_handler,
1,
inc_object,
constraints,
fe.component_mask(x_displacement) |
fe.component_mask(y_displacement) |
fe.component_mask(z_displacement));
constraints.close();
// Line below works as expected (i.e., the message "overridden value function
// called") is displayed.
//std::cout << inc_object.value(Point<dim>(0.5,0.5,0.5)) << std::endl;
return 0;
}
