132 * ********************************************************
141 *
void declare_parameters();
156 *
prm.enter_subsection (
"Global Parameters");
159 *
"Penalization parameter");
161 *
"Whether the exact solution is known");
163 *
prm.leave_subsection ();
165 *
prm.enter_subsection (
"Mesh & Refinement Parameters");
168 *
"Number identifying the domain in which we solve the problem");
170 *
"Number of global mesh refinement steps applied to initial coarse grid");
172 *
"Number of global adaptive mesh refinements");
174 *
"refinement threshold");
176 *
"coarsening threshold");
178 *
prm.leave_subsection ();
181 *
prm.enter_subsection (
"Algorithm Parameters");
184 *
"0: Preconditioned descent, 1: Newton Method");
190 *
"Maximum Number of CG iterations");
192 *
"Tolerance for CG iterations");
194 *
"Maximum Number of LS iterations");
196 *
"line search tolerance constant (c1 in Nocedal-Wright)");
198 *
"initial step length in line-search");
200 *
"Maximum Number of inner iterations");
202 *
"Threshold on norm of the derivative to declare optimality achieved");
204 *
"Threshold on norm of the derivative to declare optimality achieved in highly refined mesh");
206 *
"Number of adaptive refinement before change convergence threshold");
208 *
prm.leave_subsection ();
213 *
declare_parameters();
219 * ******************************************************************************************
231 *
virtual double value (
const Point<dim> &
pto,
const unsigned int component = 0)
const override;
262 * ******************************************************************************************
278 *
const std::vector< double > &,
286 *
virtual std::vector<std::string> get_names ()
const override;
289 *
std::vector<DataComponentInterpretation::DataComponentInterpretation>
290 *
get_data_component_interpretation ()
const override;
292 *
virtual UpdateFlags get_needed_update_flags ()
const override;
312 *
for (
unsigned int q=0;
q<n_quadrature_points; ++
q)
326 *
std::vector<std::string>
346 *
std::vector<DataComponentInterpretation::DataComponentInterpretation>
349 *
std::vector<DataComponentInterpretation::DataComponentInterpretation>
372 * ***************************************************************************************
390 *
double phi (
const double alpha)
const;
400 *
double W (
double Du2)
const;
401 *
double Wp (
double Du2)
const;
402 *
double G (
double Du2)
const;
432 *
std::string elements;
453 *
const unsigned int component = 0)
const override;
459 *
const unsigned int )
const
502 *
const unsigned int component = 0)
const override;
507 *
const unsigned int )
const
515 *
double return_value = 4.0;
516 *
return return_value;
542 *
prm.enter_subsection (
"Global Parameters");
543 *
p=prm.get_double(
"p");
544 *
prm.leave_subsection ();
545 *
prm.enter_subsection (
"Algorithm Parameters");
547 *
dir_id=prm.get_integer(
"Descent_direction");
548 *
prm.leave_subsection ();
551 *
else elements=
"P2";
564 *
dof_handler.
clear ();
582 *
std::cout <<
"iteration="<< counter+1 <<
" J(u_h)= "<<
phi_zero <<
", H1 error: "
584 *
<<
", |J'(u_h)|= "<<
system_rhs.l2_norm()<<std::endl;
588 *
std::cout <<
"iteration= " << counter+1 <<
" J(u_h)= "
644 * fill-in entry
for the solution
702 *
const unsigned int n_dofs=dof_handler.n_dofs();
739 *
dof_handler.distribute_dofs (fe);
741 *
grad_norm.reinit (dof_handler.n_dofs());
742 *
lambda.reinit (dof_handler.n_dofs());
765 *
system_matrix.reinit (sparsity_pattern);
811 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
818 *
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
822 *
cell = dof_handler.begin_active(),
823 *
endc = dof_handler.end();
824 *
for (; cell!=
endc; ++cell)
829 *
fe_values.reinit (cell);
835 *
long double coeff=0.0;
837 *
long double exponent=(p-2.0)/2*
std::log(a);
839 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
841 *
for (
unsigned int j=0;
j<dofs_per_cell; ++
j)
860 *
*fe_values.shape_value(i,
q_point)
866 *
cell->get_dof_indices (local_dof_indices);
867 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
869 *
for (
unsigned int j=0;
j<dofs_per_cell; ++
j)
870 *
system_matrix.add (local_dof_indices[i],
871 *
local_dof_indices[
j],
907 *
using FunctionMap = std::map<types::boundary_id, const Function<dim> *>;
916 *
prm.enter_subsection (
"Mesh & Refinement Parameters");
917 *
const double top_fraction=prm.get_double(
"top_fraction_of_cells");
918 *
const double bottom_fraction=prm.get_double(
"bottom_fraction_of_cells");
919 *
prm.leave_subsection ();
928 *
dof_handler.distribute_dofs(fe);
965 *
data_out.attach_dof_handler (dof_handler);
968 *
data_out.build_patches ();
969 *
std::ostringstream
p_str;
970 *
p_str << p<<
"-cycle-"<<counter;
972 *
const std::string
filename =
"solution-" +
str+
".vtu";
973 *
std::ofstream output (
filename.c_str());
974 *
data_out.write_vtu (output);
991 *
for (std::map<types::global_dof_index, double>::const_iterator
1004 *
template <
int dim>
1020 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
1024 *
std::vector<Tensor<1, dim> >
gradients(n_q_points);
1025 *
std::vector<double>
values(n_q_points);
1028 *
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1031 *
cell = dof_handler.begin_active(),
1032 *
endc = dof_handler.end();
1033 *
for (; cell!=
endc; ++cell)
1036 *
fe_values.reinit (cell);
1076 *
template <
int dim>
1087 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
1091 *
std::vector<Tensor<1, dim> >
gradients(n_q_points);
1093 *
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1096 *
cell = dof_handler.begin_active(),
1097 *
endc = dof_handler.end();
1098 *
for (; cell!=
endc; ++cell)
1101 *
fe_values.reinit (cell);
1107 *
if (
coeff<1.0e-15)
1135 *
template <
int dim>
1146 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
1150 *
std::vector<Tensor<1, dim> >
gradients(n_q_points);
1152 *
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1155 *
cell = dof_handler.begin_active(),
1156 *
endc = dof_handler.end();
1157 *
for (; cell!=
endc; ++cell)
1160 *
fe_values.reinit (cell);
1193 *
template <
int dim>
1213 *
template <
int dim>
1216 *
unsigned int it=0;
1219 *
prm.enter_subsection (
"Algorithm Parameters");
1220 *
const unsigned int max_LS_it=prm.get_integer(
"max_LS_it");
1221 *
double init_SL=prm.get_double(
"init_step_length");
1222 *
prm.leave_subsection ();
1230 *
std::cerr <<
"Step length=" <<
alpha <<
", Value= " <<
phi_alpha;
1238 *
std::cerr <<
" accepted" << std::endl;
1240 *
std::cerr <<
" rejected" ;
1259 *
std::cerr <<
", finished with steplength= "<<
alpha<<
", fcn value= "<<
phi_alpha<<std::endl;
1264 *
std::cerr <<
", max. no. of iterations reached with steplength= "<<
alpha
1265 *
<<
", fcn value= "<<
phi_alpha<<std::endl;
1285 *
template <
int dim>
1293 *
prm.enter_subsection (
"Mesh & Refinement Parameters");
1294 *
const int domain_id=prm.get_integer(
"Code for the domain");
1295 *
const int init_ref=prm.get_integer(
"No of initial refinements");
1296 *
prm.leave_subsection ();
1359 * ElastoplasticTorsion::solve(
inner_it)
1366 *
template <
int dim>
1369 *
prm.enter_subsection (
"Algorithm Parameters");
1370 *
const unsigned int max_CG_it=prm.get_integer(
"Max_CG_it");
1371 *
const double CG_tol=prm.get_double(
"CG_tol");
1372 *
prm.leave_subsection ();
1378 *
preconditioner.initialize(system_matrix,0.25);
1399 *
std::cout <<
"Not a descent direction!" <<std::endl;
1417 *
else return false;
1426 * ElastoplasticTorsion::run
1429 *
template <
int dim>
1438 *
prm.enter_subsection (
"Mesh & Refinement Parameters");
1439 *
const int adapt_ref=prm.get_integer(
"No of adaptive refinements");
1440 *
prm.leave_subsection ();
1441 *
prm.enter_subsection (
"Algorithm Parameters");
1442 *
const int max_inner=prm.get_integer(
"Max_inner");
1443 *
const double eps=prm.get_double(
"eps");
1444 *
const double hi_eps=prm.get_double(
"hi_eps");
1445 *
const int hi_th=prm.get_integer(
"hi_th");
1446 *
const double init_p=prm.get_double(
"init_p");
1447 *
const double delta_p=prm.get_double(
"delta_p");
1448 *
prm.leave_subsection ();
1449 *
prm.enter_subsection (
"Global Parameters");
1451 *
double actual_p=prm.get_double(
"p");
1452 *
prm.leave_subsection ();
1501 *
std::cout <<
"--Preparing initial condition with p="<<p<<
" iter.= " <<
global_it<<
" .-- "<< std::endl;
1515 *
ptime=timer.cpu_time();
1547 *
std::cout <<
"============ Solving problem with p=" <<p <<
" ==================" << std::endl;
1574 *
ptime=timer.cpu_time();
1586 *
std::cout <<
"finished with H1 error: " <<
H1_error <<
", dual error (L1): "
1602 *
std::cout <<
"******** Refined mesh " << cycle <<
" ********" << std::endl;
1635 *
using namespace dealii;
1636 *
using namespace nsp;
1641 *
param.read_parameters(
"EPT.prm");
1645 *
catch (std::exception &exc)
1647 *
std::cerr << std::endl << std::endl
1648 *
<<
"----------------------------------------------------" << std::endl;
1649 *
std::cerr <<
"Exception on processing: " << std::endl
1650 *
<< exc.what() << std::endl
1651 *
<<
"Aborting!" << std::endl
1652 *
<<
"----------------------------------------------------"
1659 *
std::cerr << std::endl << std::endl
1660 *
<<
"----------------------------------------------------"
1662 *
std::cerr <<
"Unknown exception!" << std::endl
1663 *
<<
"Aborting!" << std::endl
1664 *
<<
"----------------------------------------------------"
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)
unsigned int depth_console(const unsigned int n)
numbers::NumberTraits< Number >::real_type norm() const
Number linfty_norm(const Tensor< 2, dim, Number > &t)
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
#define DEAL_II_VERSION_GTE(major, minor, subminor)
typename ActiveSelector::active_cell_iterator active_cell_iterator
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.
void hyper_rectangle(Triangulation< dim, spacedim > &tria, const Point< dim > &p1, const Point< dim > &p2, const 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)
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())
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.)
void cell_residual(Vector< double > &result, const FEValuesBase< dim > &fe, const std::vector< Tensor< 1, dim > > &input, const ArrayView< const std::vector< double > > &velocity, double factor=1.)
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
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)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation