Good morning! Dr. Bangerth,
Thanks so much for your reply.
Yes, your comment exactly pointed out the problem in my declaration that
"std::vector<Point<dim>> x(dofs_per_cell)" should be created with a size
parameter.
Now I can output "x[i]" to show vertices of each cell. This is amazing!
May I ask a follow-up question?
In looping vertices of each cell, I added this:
cell_x(i) += fe_values.shape_value(i, q_index) * x(i);
however, an error shows: no match for 'operator+='(operand types are
'double' and 'dealii::Point<1, double>').
It looks to me like, "std::vector<Point<dim>> x(dofs_per_cell)" does not
want to give "Vector<double> cell_x"; but if I created as "Vector<double>
x(dofs_per_cell)", I got an error like "cannot convert 'dealii::Point<1,
double>' to 'double' in assignment".
My .cc file for 1D Laplace Eqn. is attached ready to compile, which was
made with reference to step-3. The error stops running at Line 132.
Thank you so much!
Best regards,
Judy
On Tuesday, December 7, 2021 at 9:54:19 AM UTC-5 Wolfgang Bangerth wrote:
> On 12/7/21 03:30, Judy Lee wrote:
> >
> > But how to use "Point<dim> x" for computing with other variables
> > together? Like, I cannot implement "x * fe_values.shape_value(i,
> > q_index)". However, it does not work thru if I tried like this:
> >
> > std::vector<Point<1>> x;
> > for (const unsigned int i: fe_values.dof_indices())) {
> > x[i] = cell->vertex(i);
> > }
> > std::cout << x[i] << std::endl;
> >
> > It gives:
> > make[3]: *** [CMakeFiles/run.dir/build.make:77: CMakeFiles/run]
> > Segmentation fault
> > make[2]: *** [CMakeFiles/MakeFile2:237: CMakeFiles/run.dir/all] Error 2
> > make[1]: *** [CMakeFiles/MakeFile2:244: CMakeFiles/run.dir/rule] Error 2
> > make: *** [Makefile:202: run] Error 2
>
> Judy -- can you explain *what* it is that you want to do? It is often
> easier to answer *what* questions, but you ask a *how* question.
>
> Separately, the error you show is a segmentation fault. I don't think it
> has anything to do with the line
> x[i] = cell->vertex(i)
> but instead with the fact that your vector 'x' has size zero. You need
> to set it to the correct size before you access elements 'x[i]'.
>
> Best
> W.
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: bang...@colostate.edu
> 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 dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/df049f31-9a1d-418b-b8a8-d49ab3526f4fn%40googlegroups.com.
/* ---------------------------------------------------------------------
*
* Copyright (C) 1999 - 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.
*
* ---------------------------------------------------------------------
*
* Authors: Wolfgang Bangerth, 1999,
* Guido Kanschat, 2011
*/
#include <deal.II/grid/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/grid/grid_generator.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>
using namespace dealii;
class Step3
{
public:
Step3();
void run();
private:
void make_grid();
void setup_system();
void assemble_system();
void solve();
void output_results() const;
Triangulation<1> triangulation; //edited, reset <dim = 1>;
FE_Q<1> fe; //edited, reset <dim = 1>;
DoFHandler<1> dof_handler; //edited, reset <dim = 1>;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};
Step3::Step3()
: fe(1) //Q1 element;
, dof_handler(triangulation)
{}
void Step3::make_grid()
{
//GridGenerator::hyper_cube(triangulation, -1, 1);
//triangulation.refine_global(5);
//added lines for making 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);
std::cout << "Number of active cells: " << triangulation.n_active_cells()
<< std::endl;
}
void Step3::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 Step3::assemble_system()
{
QGauss<1> quadrature_formula(fe.degree + 1); //edited, reset <dim = 1>;
FEValues<1> fe_values(fe,
quadrature_formula,
update_values | update_gradients | update_JxW_values);
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
Vector<double> cell_x(dofs_per_cell); //added lines to count in a linear function of x, "f(x) = x" to be PDE's "right-hand-side" quantity;
std::vector<Point<1>> x(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
cell_matrix = 0;
cell_rhs = 0;
cell_x = 0; //added line to count in a linear function of x, "f(x) = x" to be PDE's "right-hand-side" quantity;
//created in the same dimension as "cell_rhs";
cell->index();
for (const unsigned int q_index : fe_values.quadrature_point_indices())
{
for (const unsigned int i : fe_values.dof_indices())
for (const unsigned int j : fe_values.dof_indices())
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
for (const unsigned int i : fe_values.dof_indices())
cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
1. * // f(x_q)
fe_values.JxW(q_index)); // dx
//added lines to output a set of parameters, to find out quadrature points in local(natural) coordinate:
for (const unsigned int i : fe_values.dof_indices()){
x[i] = cell->vertex(i);
std::cout << " Print x(local_dof_indices[i]) for cell_x " << std::endl;
std::cout << Point<1>(local_dof_indices[i]) << std::endl;
std::cout << x[i] << std::endl;
cell_x(i) += fe_values.shape_value(i, q_index) * x[i]; // phi_i(x_q) * x_i; cannot run thru?...
//cell_x(i) += fe_values.shape_value(i, q_index) * 1.; // phi_i(x_q); Run thru!
}
}
cell->get_dof_indices(local_dof_indices);
for (const unsigned int i : fe_values.dof_indices())
for (const unsigned int j : fe_values.dof_indices())
system_matrix.add(local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i, j));
for (const unsigned int i : fe_values.dof_indices())
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
std::map<types::global_dof_index, double> boundary_values;
VectorTools::interpolate_boundary_values(dof_handler,
0,
Functions::ZeroFunction<1>(), //edited, reset <dim = 1>;
boundary_values);
//added lines for prescribing a Dirichlet BC @ x = L;
double g2 = 0.001; //EDIT; Functions::ConstantFunction<dim>(g2),
VectorTools::interpolate_boundary_values(dof_handler,
1,
Functions::ConstantFunction<1>(g2), //edited, reset <dim = 1>;
boundary_values);
MatrixTools::apply_boundary_values(boundary_values,
system_matrix,
solution,
system_rhs);
}
void Step3::solve()
{
SolverControl solver_control(1000, 1e-12);
SolverCG<Vector<double>> solver(solver_control);
solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
}
void Step3::output_results() const
{
DataOut<1> data_out; //edited, reset <dim = 1>;
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 Step3::run()
{
make_grid();
setup_system();
assemble_system();
solve();
output_results();
}
int main()
{
deallog.depth_console(1);
Step3 laplace_problem;
laplace_problem.run();
return 0;
}
