Dr. Bangerth,
Thank you so much for the reply that points out a quadrature point may have
one or more coordinate.
(1) I reset <dim = 1> thru step-8 and step-3. Was that sufficient to make
the program ready for 1-dimensional(1D) PDE?
If not sufficient, I might get more coordinate returned for each quadrature
point; for example, q_index = {0, 1} can be 2D coordinate for each
quadrature point, that caused error?
(2) Regarding to "obtain quadrature point coordinate locally", my .cc file
for 1D PDE with pure Dirichlet BCs is attached ready to compile for
testing, which was made with reference to step-8 that looped quadrature
points inside each cell locally:
void::assemble_system(){
...
for (const unsigned int q_index : fe_values.quadrature_point_indices())
{
std::cout << fe.degree << std::endl; //added line to output;
std::cout << q_index << std::endl; //added line to output;
std::cout << fe_values.quadrature_point << std::endl; //added line to
output;
std::cout << fe_values.JxW(q_index) << std::endl; //added line to output;
...
}
...
The added lines returned data for each cell, like:
(cell 0:)
2
0
0.00211325
0.005
2
1
0.00788675
0.005
(cell 1:)
2
0
0.0121132
0.005
2
1
0.0178860
0.005
...
The returned data from fe_values.quadrature_point {0.00211325, 0.00788675,
0.0121132, 0.0178860, ...}, look like quadrature points in physical
coordinate, moving along 1D domain of (0, 0.1) with looping cells. However,
I did not figure out more information on them.
When I switch to Q2 element, fe.degree became 2, but fe_values.JxW(q_index)
does not return expected weights of {1./6., 2./3., 1./6.}. Was that due to
Jacobians of Q2 element? I must be wrong in guessing and obtaining
quadrature point coordinate either locally or globally, as well as their
quadrature weights.
Thank you so much!
Best regards,
Judy
On Wednesday, December 1, 2021 at 8:01:22 PM UTC-5 Wolfgang Bangerth wrote:
>
> Judy,
>
> > How to obtain the value of each quadrature point locally with reference
> > to q_index, (I made this: std::cout << fe_values.quadrature_point <<, in
> > step-3 and step-8, since q_index is looping there), when loops each
> > cell, if my 1D PDE's right-hand-side is a function of x, f(x)?
> >
> > In step-3 and step-8, I set to Q1 element by switching a parameter in
> > Class Constructor. When I std::cout << q_index <<, I saw 2 indices {0,
> > 1} returned, does that mean 2 quadrature points are used? I also
> > std::cout << fe_values.JxW(q_index) <<, that gave me 0.005 corresponding
> > to both indices, so I guess 0.5 for quadrature weights, since I meshed
> > 10 cells in 1D domain of (0, 0.1). However, when I std::cout <<
> > fe_values.quadrature_point <<, I did not see expected values {0, 1}
> > locally for weights {0.5, 0.5}, when loops each cell.
> >
> > I must be wrong in guessing fe_values.JxW(q_index) and
> > fe_values.quadrature_point. Does step-3 or step-8 use 2 quadrature
> > points for Q1 element for 1D PDE, and how to obtain the value of each
> > quadrature point locally with reference to q_index, if I want to build
> > f(x), in terms of quadrature point locally?
>
> I have to admit that I don't understand your questions, principally
> because I don't understand the (pseudo)code you show. When you say
>
> (I made this: std::cout << fe_values.quadrature_point << ...
>
> what specifically does that look like? This code does not compile, so I
> don't know what it is you tried and I can't match that with the output
> you show.
>
> Can you re-state your question where you expand all of the little code
> snippets into the exact code you use?
>
> Separately: What do you mean by "quadrature point value"? A quadrature
> point is a *point*: it has one or more coordinates. The *value* is
> something that *a function returns* to you when evaluated at a point.
>
> Best
> W.
>
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: [email protected]
> www: http://www.math.colostate.edu/~bangerth/
>
--
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/8c0ba7d2-e1ef-414f-8898-b8ae7d7f6e2fn%40googlegroups.com.
/* ---------------------------------------------------------------------
*
* Copyright (C) 2000 - 2021 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE.md at
* the top level directory of deal.II.
*
* ---------------------------------------------------------------------
*
* Author: Wolfgang Bangerth, University of Heidelberg, 2000
*/
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/tensor.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/lac/affine_constraints.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/dofs/dof_handler.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>
namespace Step8
{
using namespace dealii;
template <int dim>
class ElasticProblem
{
public:
ElasticProblem();
void run();
private:
void setup_system();
void assemble_system();
void solve();
void refine_grid();
void output_results(const unsigned int cycle) const;
Triangulation<dim> triangulation;
DoFHandler<dim> dof_handler;
FESystem<dim> fe;
AffineConstraints<double> constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};
template <int dim>
void right_hand_side(const std::vector<Point<dim>> &points,
std::vector<Tensor<1, dim>> & values)
{
AssertDimension(values.size(), points.size());
Assert(dim >= 2, ExcNotImplemented());
Point<dim> point_1, point_2;
point_1(0) = 0.5;
point_2(0) = -0.5;
for (unsigned int point_n = 0; point_n < points.size(); ++point_n)
{
if (((points[point_n] - point_1).norm_square() < 0.2 * 0.2) ||
((points[point_n] - point_2).norm_square() < 0.2 * 0.2))
values[point_n][0] = 1.0;
else
values[point_n][0] = 0.0;
if (points[point_n].norm_square() < 0.2 * 0.2)
values[point_n][1] = 1.0;
else
values[point_n][1] = 0.0;
}
}
template <int dim>
ElasticProblem<dim>::ElasticProblem()
: dof_handler(triangulation)
, fe(FE_Q<dim>(1), dim)
{}
template <int dim>
void ElasticProblem<dim>::setup_system()
{
dof_handler.distribute_dofs(fe);
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
VectorTools::interpolate_boundary_values(dof_handler,
0,
Functions::ZeroFunction<dim>(dim),
constraints);
//added lines for prescribing a Dirichlet BC @ x = L: u(0) = 0, u(L) = g2;
double g2 = 0.001;
VectorTools::interpolate_boundary_values(dof_handler,
1,
Functions::ConstantFunction<dim>(g2),
constraints);
constraints.close();
DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler,
dsp,
constraints,
/*keep_constrained_dofs = */ false);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
}
template <int dim>
void ElasticProblem<dim>::assemble_system()
{
QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
const unsigned int n_q_points = quadrature_formula.size();
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);
std::vector<double> lambda_values(n_q_points);
std::vector<double> mu_values(n_q_points);
Functions::ConstantFunction<dim> lambda(1.), mu(1.);
//std::vector<Tensor<1, dim>> rhs_values(n_q_points); //edited; reset "rhs_values" to Constant at the moment, as shown in Line 192-193;
for (const auto &cell : dof_handler.active_cell_iterators())
{
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit(cell);
lambda.value_list(fe_values.get_quadrature_points(), lambda_values);
mu.value_list(fe_values.get_quadrature_points(), mu_values);
//right_hand_side(fe_values.get_quadrature_points(), rhs_values); //edited; reset "rhs_values" to Constant at the moment, as shown in Line 192-193;
for (const unsigned int i : fe_values.dof_indices())
{
const unsigned int component_i =
fe.system_to_component_index(i).first;
for (const unsigned int j : fe_values.dof_indices())
{
const unsigned int component_j =
fe.system_to_component_index(j).first;
for (const unsigned int q_point :
fe_values.quadrature_point_indices())
{
//added lines to output a set of values, to find out quadrature points in local(natural) or physical coordinate:
std::cout << fe.degree + 1 << std::endl; //with reference to Line 117;
std::cout << q_point << std::endl;
std::cout << fe_values.quadrature_point(q_point) << std::endl;
std::cout << fe_values.JxW(q_point) << std::endl;
cell_matrix(i, j) +=
(
(fe_values.shape_grad(i, q_point)[component_i] *
fe_values.shape_grad(j, q_point)[component_j] *
lambda_values[q_point])
+
(fe_values.shape_grad(i, q_point)[component_j] *
fe_values.shape_grad(j, q_point)[component_i] *
mu_values[q_point])
+
((component_i == component_j) ?
(fe_values.shape_grad(i, q_point) *
fe_values.shape_grad(j, q_point) *
mu_values[q_point]) :
0)
) *
fe_values.JxW(q_point);
}
}
}
for (const unsigned int i : fe_values.dof_indices())
{
const unsigned int component_i =
fe.system_to_component_index(i).first;
for (const unsigned int q_point :
fe_values.quadrature_point_indices())
cell_rhs(i) += fe_values.shape_value(i, q_point) *
//rhs_values[q_point][component_i] *
1.0 * //edited; reset "rhs_values" to Constant at the moment;
fe_values.JxW(q_point);
}
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(
cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
}
}
template <int dim>
void ElasticProblem<dim>::solve()
{
SolverControl solver_control(1000, 1e-12);
SolverCG<Vector<double>> cg(solver_control);
PreconditionSSOR<SparseMatrix<double>> preconditioner;
preconditioner.initialize(system_matrix, 1.2);
cg.solve(system_matrix, solution, system_rhs, preconditioner);
constraints.distribute(solution);
}
template <int dim>
void ElasticProblem<dim>::refine_grid()
{
Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
KellyErrorEstimator<dim>::estimate(dof_handler,
QGauss<dim - 1>(fe.degree + 1),
{},
solution,
estimated_error_per_cell);
GridRefinement::refine_and_coarsen_fixed_number(triangulation,
estimated_error_per_cell,
0.3,
0.03);
triangulation.execute_coarsening_and_refinement();
}
template <int dim>
void ElasticProblem<dim>::output_results(const unsigned int cycle) const
{
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
std::vector<std::string> solution_names;
switch (dim)
{
case 1:
solution_names.emplace_back("displacement");
break;
case 2:
solution_names.emplace_back("x_displacement");
solution_names.emplace_back("y_displacement");
break;
case 3:
solution_names.emplace_back("x_displacement");
solution_names.emplace_back("y_displacement");
solution_names.emplace_back("z_displacement");
break;
default:
Assert(false, ExcNotImplemented());
}
data_out.add_data_vector(solution, solution_names);
data_out.build_patches();
std::ofstream output("solution-" + std::to_string(cycle) + ".vtk");
data_out.write_vtk(output);
}
template <int dim>
void ElasticProblem<dim>::run()
{
for (unsigned int cycle = 0; cycle < 1; ++cycle) //edited, reset <"cycle < 1">;
{
std::cout << "Cycle " << cycle << ':' << std::endl;
if (cycle == 0)
{
//GridGenerator::hyper_cube(triangulation, -1, 1);
//triangulation.refine_global(4);
//added lines to make 1D Physical Domain (0, 0.1);
double L = 0.1;
double x_min = 0.;
double x_max = L;
const unsigned int numberOfCells = 10;
const unsigned int meshDimensions = numberOfCells;
GridGenerator::subdivided_hyper_cube(triangulation, meshDimensions, x_min, x_max);
}
else
refine_grid();
std::cout << " Number of active cells: "
<< triangulation.n_active_cells() << std::endl;
setup_system();
std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
assemble_system();
solve();
output_results(cycle);
}
}
} // namespace Step8
int main()
{
try
{
Step8::ElasticProblem<1> elastic_problem_2d; //edited, reset <dim = 1>;
elastic_problem_2d.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
