Hello Jean,
First of all, thanks a lot for taking time to write to me and explaining
this issue.
I am also very keen on solving this issue.
I looked at the code you sent and indeed it looks like what you say is
correct.
However, If I calculate the gradients, it is giving wrong results and
giving gradients on the mesh before movement.
I took the gradients before and after moving the mesh and it is giving me
the same values which has been happening with me.
I am also attaching my file and CMakelists.txt for you to review.
I am running in debug mode but the issue is still unresolved.
I am taking the gradient of "f" which I initialized to the material
coordinates.
Is this a bug in the library or the version that I am using (8.3.2) ?
gradient at qpoint 0 before mesh movement:: 1 0 0 1
Volume before mesh movement: Using GridTools: 1 Calculated: 1
gradient at qpoint 0 after mesh movement :: 1 0 0 1
Volume after mesh movement: Using GridTools: 0.95 Calculated: 1
gradient at qpoint 0 after mesh movement and with new fe_values object ::
1.05263 0 0 1
On Tuesday, May 31, 2016 at 3:19:13 PM UTC-4, Jean-Paul Pelteret wrote:
>
> Hi Rajat,
>
> Without your code, its difficult to say what's going wrong here. It should
> (and does) work as I've explained. In summary, FEValues gets to know about
> any changes you've made to the underlying because its reinit() call takes
> an *active* cell which always knows its geometry (you manipulate this
> directly when you move the mesh).
>
> I've written a small test case to prove to myself that your "code 1"
> should work as expected, and indeed this is the case. I've attached it for
> you to investigate, and this is the expected result:
> Volume before mesh movement: Using GridTools: 1 Calculated: 1
> Volume after mesh movement: Using GridTools: 0.75 Calculated: 0.75
>
> One thing that *might* be the problem is that you might be iterating over
> inactive cells with "dof_handler.begin()" instead of
> "dof_handler.begin_active()". If you only have a single cell and
> accidentally do this, then even in debug mode you do get a result but it is
> incorrect:
> Volume before mesh movement: Using GridTools: 1 Calculated: 1
> Volume after mesh movement: Using GridTools: 0.75 Calculated: 1
>
> However, once you perform some refinement then this error would be picked
> up in debug mode (are you running your tests in debug mode?).
>
> Does that help? If you solve the issue, then I'd be interested to know
> what the problem was.
> J-P
>
> On Tuesday, May 31, 2016 at 4:17:56 PM UTC+2, RAJAT ARORA wrote:
>>
>> Hello Jean,
>>
>> Thanks for your reply.
>>
>> In my case, since I am physically moving the mesh ( as I was discussing
>> here <https://groups.google.com/forum/#!topic/dealii/8I21UqIWAmU>), even
>> after calling
>> fe_values.reinit(cell), I am getting the gradients on the geometry before
>> mesh movement which should not be the case.
>>
>>
>>
>> Also, note that once I move the mesh, this happens only if I use the same
>> fe_values object. But if I create a new FEValues object (lets call it
>> fe_values_new), then if I use fe_values_new.get_function_gradients(), I get
>> the gradients on the updated geometry. Code given below for reference.
>>
>>
>> Code 1: (Does not work)
>> // Everything is in one big function
>>
>> FEValues<dim> fe_values (fe, quadrature_formula,
>> update_values
>> <https://dealii.org/8.4.1/doxygen/deal.II/group__feaccess.html#ggaa94b67d2fdcc390690c523f28019e52fa4057ca2f127aa619c65886a9d3ad4aea>
>>
>> | update_gradients
>> <https://dealii.org/8.4.1/doxygen/deal.II/group__feaccess.html#ggaa94b67d2fdcc390690c523f28019e52facbcc430975fa6af05f75ca786dc6fe20>
>>
>> |
>> update_quadrature_points
>> <https://dealii.org/8.4.1/doxygen/deal.II/group__feaccess.html#ggaa94b67d2fdcc390690c523f28019e52fad5c9ff886b9615349a5d04a6c782df4a>
>>
>> | update_JxW_values
>> <https://dealii.org/8.4.1/doxygen/deal.II/group__feaccess.html#ggaa94b67d2fdcc390690c523f28019e52fa714204722e9eeb43aadbd0d5ddc48c85>
>> );
>>
>> // do some stuff
>>
>> // move mesh like step 17, adding displacements to cell->vertex(i)
>>
>> fe_values.reinit(cell);
>>
>> fe_values.get_function_gradients();
>>
>>
>>
>> Code 2:(It works)
>> // Everything is in one big function
>>
>> FEValues<dim> fe_values = new FEValues<dim>(fe, quadrature_formula,
>> update_values
>> <https://dealii.org/8.4.1/doxygen/deal.II/group__feaccess.html#ggaa94b67d2fdcc390690c523f28019e52fa4057ca2f127aa619c65886a9d3ad4aea>
>> | update_gradients
>> <https://dealii.org/8.4.1/doxygen/deal.II/group__feaccess.html#ggaa94b67d2fdcc390690c523f28019e52facbcc430975fa6af05f75ca786dc6fe20>
>> |
>> update_quadrature_points
>> <https://dealii.org/8.4.1/doxygen/deal.II/group__feaccess.html#ggaa94b67d2fdcc390690c523f28019e52fad5c9ff886b9615349a5d04a6c782df4a>
>> | update_JxW_values
>> <https://dealii.org/8.4.1/doxygen/deal.II/group__feaccess.html#ggaa94b67d2fdcc390690c523f28019e52fa714204722e9eeb43aadbd0d5ddc48c85>
>> );
>>
>> // do some stuff
>>
>> // move mesh like step 17, adding displacements to cell->vertex(i)
>>
>> delete fe_values;
>>
>> fe_values = new FEValues<dim> (fe, quadrature_formula,
>> update_values
>> <https://dealii.org/8.4.1/doxygen/deal.II/group__feaccess.html#ggaa94b67d2fdcc390690c523f28019e52fa4057ca2f127aa619c65886a9d3ad4aea>
>> | update_gradients
>> <https://dealii.org/8.4.1/doxygen/deal.II/group__feaccess.html#ggaa94b67d2fdcc390690c523f28019e52facbcc430975fa6af05f75ca786dc6fe20>
>> |
>> update_quadrature_points
>> <https://dealii.org/8.4.1/doxygen/deal.II/group__feaccess.html#ggaa94b67d2fdcc390690c523f28019e52fad5c9ff886b9615349a5d04a6c782df4a>
>> | update_JxW_values
>> <https://dealii.org/8.4.1/doxygen/deal.II/group__feaccess.html#ggaa94b67d2fdcc390690c523f28019e52fa714204722e9eeb43aadbd0d5ddc48c85>
>> );
>>
>> fe_values.reinit(cell);
>>
>> fe_values.get_function_gradients();
>>
>>
>>
>>
>>
>>
>> Does that mean the instant FEValues object is created, it takes that as a
>> reference configuration and do all the gradients on this reference
>> configuration ?
>>
>> I also read the posts related to this issue here
>> <https://groups.google.com/forum/#!searchin/dealii/mesh$20fevalues/dealii/KvqA-fQAYgM/yQdzNZsRSOUJ>
>> and here
>> <https://groups.google.com/forum/#!searchin/dealii/mesh$20fevalues/dealii/B2v7qSLrnck/0e44LC4bCQAJ>
>> .
>>
>>
>>
>>
>>
>> On Tuesday, May 31, 2016 at 1:41:10 AM UTC-4, Jean-Paul Pelteret wrote:
>>>
>>> If you physically move the mesh vertices, like you were discussing here
>>> <https://groups.google.com/forum/#!topic/dealii/8I21UqIWAmU>, then when
>>> you call fe_values.reinit(cell), the mapping between the current cell
>>> geometry (as defined by the triangulation) and isoparametric cell are
>>> accounted for and used when computing gradient quantities.
>>>
>>> For example, compare the update_quadrature_point_history functions in
>>> step-18
>>> <https://dealii.org/8.4.0/doxygen/deal.II/step_18.html#TopLevelupdate_quadrature_point_history>
>>>
>>> and step-44
>>> <https://dealii.org/8.4.0/doxygen/deal.II/step_44.html#Solidupdate_qph_incremental>.
>>>
>>> In step-44 the mesh is not physically moved, so "get_function_gradients"
>>> applied to the total displacement vector will return elements of
>>> \frac{\partial \mathbf{u}}{\partial \mathbf{X}} (i.e. the material
>>> gradient), while in step-18 the mesh is moved and "get_function_gradients"
>>> returns elements of \frac{\partial \mathbf{u}}{\partial \mathbf{x}} (i.e.
>>> the spatial gradient). However, you will note that in step-18, only the
>>> incremental displacement gradient is computed as the stresses are
>>> systematically incremented using this result.
>>>
>>> On Tuesday, May 31, 2016 at 2:54:32 AM UTC+2, RAJAT ARORA wrote:
>>>>
>>>> Hello all,
>>>>
>>>> I have a question about FEValues<dim> object that we use in the
>>>> computation of element stiffness and other places.
>>>>
>>>>
>>>> If I declare fe_values object and get then move the mesh and then try
>>>> to call fe_values.get_function_gradients()
>>>> will this give me gradients on the deformed mesh or the mesh that was
>>>> chosen to begin with ? (Code structure below)
>>>>
>>>>
>>>> FEValues<dim> fe_values (fe, quadrature_formula,
>>>> update_values
>>>> <https://dealii.org/8.4.1/doxygen/deal.II/group__feaccess.html#ggaa94b67d2fdcc390690c523f28019e52fa4057ca2f127aa619c65886a9d3ad4aea>
>>>>
>>>> | update_gradients
>>>> <https://dealii.org/8.4.1/doxygen/deal.II/group__feaccess.html#ggaa94b67d2fdcc390690c523f28019e52facbcc430975fa6af05f75ca786dc6fe20>
>>>>
>>>> |
>>>> update_quadrature_points
>>>> <https://dealii.org/8.4.1/doxygen/deal.II/group__feaccess.html#ggaa94b67d2fdcc390690c523f28019e52fad5c9ff886b9615349a5d04a6c782df4a>
>>>>
>>>> | update_JxW_values
>>>> <https://dealii.org/8.4.1/doxygen/deal.II/group__feaccess.html#ggaa94b67d2fdcc390690c523f28019e52fa714204722e9eeb43aadbd0d5ddc48c85>
>>>> );
>>>>
>>>> // do some stuff
>>>>
>>>> // move mesh like step 17, adding displacements to cell->vertex(i)
>>>>
>>>> fe_values.get_function_gradients();
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> If it gives gradients on the initial mesh, then how does it know those
>>>> coordinates once you change the coordinates of the
>>>> vertices of the triangulation ?
>>>>
>>>>
>>>>
>>>
--
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].
For more options, visit https://groups.google.com/d/optout.
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/base/quadrature_lib.h>
#include <iostream>
#include <fstream>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/lac/vector.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_q.h>
#include <fstream>
#include <iostream>
using namespace dealii;
template<int dim>
double
integral_volume (const DoFHandler<dim> &dof_handler,
FEValues<dim> *fe_values)
{
double vol = 0.0;
for (typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active();
cell != dof_handler.end();
++cell)
{
fe_values->reinit(cell);
for (unsigned int qp_cell=0; qp_cell<fe_values->get_quadrature().size(); ++qp_cell)
vol += fe_values->JxW(qp_cell);
}
return vol;
}
template <int dim>
class velocityeqn
{
public:
velocityeqn(Triangulation<dim> &tria)
:
fe(FE_Q<dim>(1), dim),
dof_handler(tria),
quadrature_formula(dim+1)
{
dof_handler.distribute_dofs(fe);
}
void move_mesh()
{
for (typename DoFHandler<dim>::active_cell_iterator
cell = this->dof_handler.begin_active();
cell != this->dof_handler.end();
++cell)
{
// Decrease the x-length of the geometry by 25%
// By moving only the right-most vertices
// Note: This distorts only the right two
// cells
for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
if (std::fabs(cell->vertex(v)[0] - 1.0) < 1e-9)
cell->vertex(v)[0] *= 0.95;
}
}
~velocityeqn()
{
dof_handler.clear();
}
Vector<double> velocity;
FESystem<dim> fe;
DoFHandler<dim> dof_handler;
QGauss<dim> quadrature_formula;
};
template <int dim>
class Plasticeqn
{
public:
Plasticeqn(Triangulation<dim> &tria)
:
fe(FE_Q<dim>(1), dim),
dof_handler(tria),
quadrature_formula(dim+1)
{
dof_handler.distribute_dofs(fe);
f.reinit(dof_handler.n_dofs());
typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
endc = dof_handler.end();
for(;cell != endc; ++cell)
{
for(int v = 0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
for(int d = 0; d<dim; ++d)
{
// std::cout << "v " << v << "d " << d <<std::endl;
f(cell->vertex_dof_index(v,d)) = cell->vertex(v)[d];
}
}
fe_values = new FEValues<dim> (fe, quadrature_formula,
update_JxW_values | update_gradients);
}
~Plasticeqn()
{
delete fe_values;
dof_handler.clear();
}
void clear()
{
delete fe_values;
fe_values = new FEValues<dim> (fe, quadrature_formula,
update_JxW_values | update_gradients);
}
Tensor<2,dim> get_gradient() // returns gradient at 0th quadrature point
{
int n_q_points = this->quadrature_formula.size();
std::vector<std::vector< Tensor< 1, dim > > >
dfdx(n_q_points, std::vector<Tensor<1, dim> >(dim));
for (typename DoFHandler<dim>::active_cell_iterator
cell = this->dof_handler.begin_active();
cell != this->dof_handler.end();
++cell)
{
this->fe_values->reinit(cell);
this->fe_values->get_function_gradients(this->f, dfdx);
}
Tensor<2,dim> mygradf;
for (int i = 0; i < dim; ++i)
for (int j = 0; j < dim; ++j)
{
mygradf[i][j] = dfdx[0][i][j];
}
return mygradf;
}
FESystem<dim> fe;
DoFHandler<dim> dof_handler;
QGauss<dim> quadrature_formula;
Vector<double> f;
FEValues<dim> *fe_values;
};
int main()
{
const int dim = 2;
Triangulation<dim> tria;
const int degree = 1;
FE_Q<dim> fe(degree);
QGauss<dim> qf_cell(degree+1);
DoFHandler<dim> dof_handler(tria);
// Single cell of unit length
GridGenerator::hyper_cube(tria, 0.0, 1.0);
tria.refine_global(0);
velocityeqn<dim>* veqn = new velocityeqn<dim>( tria);
Plasticeqn<dim>* feqn = new Plasticeqn<dim>( tria);
dof_handler.distribute_dofs(fe);
std::cout << "gradient at qpoint 0 before mesh movement:: " << feqn->get_gradient() << std::endl;
FEValues<dim> *fe_values = new FEValues<dim> (fe, qf_cell,
update_JxW_values);
// Compute volume of the grid by integration
{
std::cout
<< "Volume before mesh movement: "
<< " Using GridTools: " << GridTools::volume(tria)
<< " Calculated: " << integral_volume(dof_handler, fe_values)
<< std::endl;
std::ofstream output_file ("mesh-initial.eps");
GridOut().write_eps(tria,output_file);
}
veqn->move_mesh();
std::cout << "gradient at qpoint 0 after mesh movement :: " << feqn->get_gradient() << std::endl;
// Move one of the grid vertices
// Now we measure the volume again.
// Since the actual geometry has been updates
// (i.e. the "reference" configuration has changed)
// the computed volume will change as well
{
std::cout
<< "Volume after mesh movement: "
<< " Using GridTools: " << GridTools::volume(tria)
<< " Calculated: " << integral_volume(dof_handler, fe_values)
<< std::endl;
std::ofstream output_file ("mesh-final.eps");
GridOut().write_eps(tria,output_file);
}
feqn->clear();
std::cout << "gradient at qpoint 0 after mesh movement and with new fe_values object :: " << feqn->get_gradient() << std::endl;
delete fe_values;
delete veqn;
delete feqn;
return 0;
}
##
# CMake script for the fdm tutorial program:
##
# Set the name of the project and target:
SET(TARGET "main")
add_compile_options(-Wno-comment -Wno-unused)
# Declare all source files the target consists of:
include_directories("${PROJECT_SOURCE_DIR}/src")
FILE(GLOB SOURCES "./src/*.cc")
SET(TARGET_SRC ./mesh_movement.cc)
# Usually, you will not need to modify anything beyond this point...
CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
FIND_PACKAGE(deal.II 8.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 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()
#
# Are all dependencies fullfilled?
#
IF(NOT DEAL_II_WITH_PETSC OR NOT DEAL_II_WITH_P4EST)
MESSAGE(FATAL_ERROR "
Error! The deal.II library found at ${DEAL_II_PATH} was not configured with
DEAL_II_WITH_PETSC = ON
DEAL_II_WITH_P4EST = ON
One or all of these are OFF in your installation but are required for this
tutorial step."
)
ENDIF()
DEAL_II_INITIALIZE_CACHED_VARIABLES()
PROJECT(${TARGET} CXX)
DEAL_II_INVOKE_AUTOPILOT()
