590 *
if (p.norm_square() < 0.25)
611 *
if (p.norm_square() < 0.25)
622 * <a name=
"step_65-ThePoissonProblemclass"></a>
659 *
Vector<double> solution;
687 * <a name=
"step_65-Gridcreationandinitializationofthemanifolds"></a>
735 * `cell->face(f)->set_all_manifold_ids(1)`
to set
the manifold
id
764 *
if (
std::
abs(face->vertex(v).norm_square() - 0.25) > 1
e-12)
768 *
face->set_all_manifold_ids(1);
770 *
if (cell->center().norm_square() < 0.25)
771 *
cell->set_material_id(1);
773 *
cell->set_material_id(0);
833 *
dof_handler.distribute_dofs(fe);
834 *
std::cout <<
" Number of active cells: "
836 *
std::cout <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
842 *
constraints.
clear();
848 *
constraints.close();
854 *
sparsity_pattern.copy_from(
dsp);
855 *
system_matrix.reinit(sparsity_pattern);
857 *
solution.reinit(dof_handler.n_dofs());
865 * <a name=
"step_65-Assemblyofthesystemmatrixandrighthandside"></a>
872 * set
the number
of quadrature points
to the polynomial degree plus
937 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
943 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
945 *
for (
const auto &cell : dof_handler.active_cell_iterators())
948 *
fe_values.reinit(cell);
950 *
for (
unsigned int q_index = 0; q_index < n_q_points; ++q_index)
953 *
coefficient(fe_values.quadrature_point(q_index));
954 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
956 *
for (
unsigned int d = 0;
d < dim; ++
d)
959 *
fe_values.shape_grad(i, q_index)[
d];
961 *
(fe_values.shape_value(i, q_index) *
963 *
fe_values.JxW(q_index));
969 *
cell->get_dof_indices(local_dof_indices);
970 *
constraints.distribute_local_to_global(
980 * <a name=
"step_65-Solutionofthelinearsystem"></a>
998 *
preconditioner.initialize(system_matrix);
1000 *
solver.solve(system_matrix, solution,
system_rhs, preconditioner);
1001 *
constraints.distribute(solution);
1003 *
std::cout <<
" Number of solver iterations: "
1004 *
<< solver_control.last_step() << std::endl;
1012 * <a name=
"step_65-Outputofthesolutionandcomputationoferrors"></a>
1045 *
href=
"https://github.com/dealii/dealii/wiki/Notes-on-visualizing-high-order-output">
wiki
1072 *
flags.write_higher_order_cells =
true;
1073 *
data_out.set_flags(flags);
1075 *
data_out.attach_dof_handler(dof_handler);
1076 *
data_out.add_data_vector(solution,
"solution");
1079 *
for (
const auto &cell :
triangulation.active_cell_iterators())
1080 *
material_ids[cell->active_cell_index()] = cell->
material_id();
1081 *
data_out.add_data_vector(material_ids,
"material_ids");
1083 *
data_out.build_patches(mapping,
1087 *
std::ofstream file(
1089 *
std::to_string(
triangulation.n_global_levels() - 10 + 2 * dim) +
1093 *
data_out.write_vtu(file);
1124 *
std::cout <<
" L2 error vs exact solution: "
1134 *
std::cout <<
" H1 error vs exact solution: "
1161 *
std::cout <<
" Max cell-wise error estimate: "
1171 * <a name=
"step_65-ThePoissonProblemrunfunction"></a>
1172 * <
h3>
The PoissonProblem::run() function</
h3>
1202 *
std::cout << std::endl
1203 *
<<
"====== Running with the basic MappingQ class ====== "
1211 *
postprocess(mapping);
1213 *
timer.print_summary();
1233 *
<<
"====== Running with the optimized MappingQCache class ====== "
1242 *
std::cout <<
" Memory consumption cache: "
1248 *
postprocess(mapping);
1250 *
timer.print_summary();
1268<a name=
"step_65-Programoutput"></a><
h3>
Program output</
h3>
1292+---------------------------------------------+------------+------------+
1296+---------------------------------+-----------+------------+------------+
1298|
Compute constraints | 1 | 0.109s | 0.22% |
1302|
Write output | 1 | 4.85s | 9.8% |
1303+---------------------------------+-----------+------------+------------+
1316+---------------------------------------------+------------+------------+
1320+---------------------------------+-----------+------------+------------+
1322|
Compute constraints | 1 | 0.00336s | 0% |
1325|
Initialize mapping cache | 1 | 4.96s | 27% |
1327|
Write output | 1 | 0.875s | 4.8% |
1328+---------------------------------+-----------+------------+------------+
1336(solution
or right
hand side) vector,
ObserverPointer< const Triangulation< dim, spacedim > > triangulation
virtual void build_patches(const unsigned int n_subdivisions=0)
void mTmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
void initialize(const Mapping< dim, spacedim > &mapping, const Triangulation< dim, spacedim > &triangulation)
Abstract base class for mapping classes.
Number linfty_norm(const Tensor< 2, dim, Number > &t)
static constexpr std::size_t memory_consumption()
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
const Triangulation< dim, spacedim > * triangulation
void initialize(const Triangulation< dim, spacedim > &triangulation)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
void hyper_shell(Triangulation< dim, spacedim > &tria, const Point< spacedim > ¢er, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, bool colorize=false)
void hyper_ball(Triangulation< dim, spacedim > &tria, const Point< spacedim > ¢er={}, const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void merge_triangulations(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result, const double duplicated_vertex_tolerance=1.0e-12, const bool copy_manifold_ids=false, const bool copy_boundary_ids=false)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
@ matrix
Contents is actually a matrix.
constexpr types::blas_int one
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T sum(const T &t, const MPI_Comm mpi_communicator)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
DEAL_II_HOST constexpr TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
int(&) functions(const void *v1, const void *v2)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation