770 *
sparsity_pattern.copy_from(
dsp);
772 *
system_matrix.reinit(sparsity_pattern);
779 * <a name=
"step_8-ElasticProblemassemble_system"></a>
813 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
819 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
830 *
std::vector<double>
mu_values(n_q_points);
850 *
std::vector<Tensor<1, dim>>
rhs_values(n_q_points);
857 *
for (
const auto &cell : dof_handler.active_cell_iterators())
871 *
mu.value_list(fe_values.get_quadrature_points(),
mu_values);
882 * <
code>fe.system_to_component_index(i).first</
code> function call
904 *
for (
const unsigned int i : fe_values.dof_indices())
907 *
fe.system_to_component_index(i).
first;
909 *
for (
const unsigned int j : fe_values.dof_indices())
912 *
fe.system_to_component_index(
j).
first;
914 *
for (
const unsigned int q_point :
915 *
fe_values.quadrature_point_indices())
959 *
(fe_values.shape_grad(i,
q_point) *
960 *
fe_values.shape_grad(
j,
q_point) *
975 *
for (
const unsigned int i : fe_values.dof_indices())
978 *
fe.system_to_component_index(i).
first;
980 *
for (
const unsigned int q_point :
981 *
fe_values.quadrature_point_indices())
995 *
cell->get_dof_indices(local_dof_indices);
996 *
constraints.distribute_local_to_global(
1006 * <a name=
"step_8-ElasticProblemsolve"></a>
1007 * <
h4>ElasticProblem::solve</
h4>
1017 *
template <
int dim>
1024 *
preconditioner.initialize(system_matrix, 1.2);
1026 *
cg.solve(system_matrix, solution,
system_rhs, preconditioner);
1028 *
constraints.distribute(solution);
1035 * <a name=
"step_8-ElasticProblemrefine_grid"></a>
1036 * <
h4>ElasticProblem::refine_grid</
h4>
1052 *
template <
int dim>
1075 * <a name=
"step_8-ElasticProblemoutput_results"></a>
1087 *
To do this,
the DataOut::add_vector() function
wants a vector
of
1124 *
data_out.attach_dof_handler(dof_handler);
1158 *
data_out.build_patches();
1160 *
std::ofstream output(
"solution-" + std::to_string(cycle) +
".vtk");
1161 *
data_out.write_vtk(output);
1169 * <a name=
"step_8-ElasticProblemrun"></a>
1170 * <
h4>ElasticProblem::run</
h4>
1203 * Triangulation::refine_and_coarsen_fixed_number()
will not flag any cells
1221 *
for (
unsigned int cycle = 0; cycle < 8; ++cycle)
1223 *
std::cout <<
"Cycle " << cycle <<
':' << std::endl;
1233 *
std::cout <<
" Number of active cells: "
1238 *
std::cout <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
1251 * <a name=
"step_8-Thecodemaincodefunction"></a>
1268 *
catch (std::exception &exc)
1270 *
std::cerr << std::endl
1272 *
<<
"----------------------------------------------------"
1274 *
std::cerr <<
"Exception on processing: " << std::endl
1275 *
<< exc.what() << std::endl
1276 *
<<
"Aborting!" << std::endl
1277 *
<<
"----------------------------------------------------"
1284 *
std::cerr << std::endl
1286 *
<<
"----------------------------------------------------"
1288 *
std::cerr <<
"Unknown exception!" << std::endl
1289 *
<<
"Aborting!" << std::endl
1290 *
<<
"----------------------------------------------------"
1310<
img src=
"https://www.dealii.org/images/steps/developer/step-8.x.png" alt=
"">
1313<
img src=
"https://www.dealii.org/images/steps/developer/step-8.y.png" alt=
"">
1325this is so, remember that we have just defined our finite element as a
1326collection of two components (in <code>dim=2</code> dimensions). Nowhere have
1327we said that this is not just a pressure and a concentration (two scalar
1328quantities) but that the two components actually are the parts of a
1329vector-valued quantity, namely the displacement. Absent this knowledge, the
1330DataOut class assumes that all individual variables we print are separate
1331scalars, and VisIt and Paraview then faithfully assume that this is indeed what it is. In
1332other words, once we have written the data as scalars, there is nothing in
1333these programs that allows us to paste these two scalar fields back together as a
1334vector field. Where we would have to attack this problem is at the root,
1335namely in <code>ElasticProblem::output_results()</code>. We won't
do so here but
1338anyway that would show how this would look if implemented as discussed in
1339@ref step_22 "step-22". The vector field then looks like this (VisIt and Paraview
1340randomly select a few
1341hundred vertices from which to draw the vectors; drawing them from each
1342individual vertex would make the picture unreadable):
1344<img src="https://www.dealii.org/images/steps/developer/step-8.vectors.png" alt="">
1347We note that one may have intuitively expected the
1348solution to be symmetric about the @f$x@f$- and @f$y@f$-axes since the @f$x@f$- and
1349@f$y@f$-forces are symmetric with respect to these axes. However, the force
1350considered as a vector is not symmetric and consequently neither is
1354<a name="step_8-PlainProg"></a>
1355<h1> The plain program</h1>
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)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
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())
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
#define DEAL_II_NOT_IMPLEMENTED()
std::vector< index_type > data
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
@ general
No special properties.
constexpr types::blas_int zero
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.)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
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)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation