245#include "allheaders.h"
246#include "nonlinear_heat.h"
270 const unsigned int dofs_per_cell = fe.dofs_per_cell;
294 std::vector<types::global_dof_index> local_dof_indices(
300 std::vector<double>
consol(n_q_points);
345 for (
const auto &cell: dof_handler.active_cell_iterators())
348 fe_values.reinit(cell);
349 cell->get_dof_indices(local_dof_indices);
350 const unsigned int n_independent_variables = local_dof_indices.size();
351 const unsigned int n_dependent_variables = dofs_per_cell;
368 fe_values[t].get_function_values_from_local_dof_values(
dof_values_ad,
370 fe_values[t].get_function_gradients_from_local_dof_values(
dof_values_ad,
386 for (
unsigned int q_index = 0; q_index < n_q_points; ++q_index)
388 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
452 cell->get_dof_indices(local_dof_indices);
454 for (
unsigned int i =0;i < dofs_per_cell; ++i)
455 residual(local_dof_indices[i])+=
cell_rhs(i);
465 std::cout <<
" The Norm is :: = " << residual.l2_norm() << std::endl;
472####
The compute_jacobian() function
477 ad_helper.compute_linearization(cell_matrix);
483### General ideas involved in solving coupled nonlinear equations using Newton Raphson's technique
491We
need to solve
the set
equations immediately above using Newton-
Raphson's technique. Now suppose we have 4 noded linear element, with values as @f$p_1^{{s+1}},p_2^{{s+1}},p_3^{{s+1}},p_4^{{s+1}}@f$ and @f$T_1^{{s+1}},T_2^{{s+1}},T_3^{{s+1}},T_4^{{s+1}}@f$, then actually @f$f_i^1@f$ and @f$f_i^2@f$ are
493 f_i^1(p_1^{{s+1}},p_2^{{s+1}},p_3^{{s+1}},p_4^{{s+1}},T_1^{{s+1}},T_2^{{s+1}},T_3^{{s+1}},T_4^{{s+1}})=0 \\
494 f_i^2(p_1^{{s+1}},p_2^{{s+1}},p_3^{{s+1}},p_4^{{s+1}},T_1^{{s+1}},T_2^{{s+1}},T_3^{{s+1}},T_4^{{s+1}})=0
507 %================================================================
548 %================================================================
578 %================================================================
627 ||(function)^
k||<\delta^0 ||function^0||
642 additional_data.max_iter = 100;
669 const double tolerance) {
670 solve(
rhs, dst, tolerance);
686 additional_data.max_iter = 100;
718 const double tolerance) {
719 solve(
rhs, dst, tolerance);
729 const double tolerance) {
730 solve(
rhs, dst, tolerance);
763### Evaluating the function at previous iterations
802<a name=
"ann-mesh/README.md"></a>
817<a name=
"ann-include/allheaders.h"></a>
924<a name=
"ann-include/nonlinear_heat.h"></a>
962 *
const double delta_t;
966 *
const double alpha;
1081<a name=
"ann-nonlinear_heat.cc"></a>
1107 *
catch (std::exception &exc)
1109 *
std::cerr << std::endl
1111 *
<<
"----------------------------------------------------"
1113 *
std::cerr <<
"Exception on processing: " << std::endl
1114 *
<< exc.what() << std::endl
1115 *
<<
"Aborting!" << std::endl
1116 *
<<
"----------------------------------------------------"
1123 *
std::cerr << std::endl
1125 *
<<
"----------------------------------------------------"
1127 *
std::cerr <<
"Unknown exception!" << std::endl
1128 *
<<
"Aborting!" << std::endl
1129 *
<<
"----------------------------------------------------"
1139<a name=
"ann-source/boundary_values.cc"></a>
1163 *
double Boundary_values_left::value(
const Point<2> & ,
const unsigned int )
const
1183<a name=
"ann-source/compute_jacobian.cc"></a>
1225 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
1236 *
system_matrix = 0.0;
1237 *
std::vector<types::global_dof_index> local_dof_indices(
1240 *
std::vector<double>
consol(n_q_points);
1249 *
for (
const auto &cell: dof_handler.active_cell_iterators())
1253 *
fe_values.reinit(cell);
1254 *
cell->get_dof_indices(local_dof_indices);
1255 *
const unsigned int n_independent_variables = local_dof_indices.size();
1256 *
const unsigned int n_dependent_variables = dofs_per_cell;
1260 *
fe_values[t].get_function_values_from_local_dof_values(
dof_values_ad,
1262 *
fe_values[t].get_function_gradients_from_local_dof_values(
dof_values_ad,
1267 *
std::vector<ADNumberType>
residual_ad(n_dependent_variables,
1269 *
for (
unsigned int q_index = 0; q_index < n_q_points; ++q_index)
1271 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1288 *
if (cell->face(face_number)->boundary_id() == 3)
1290 *
fe_face_values.reinit(cell, face_number);
1293 *
for (
unsigned int i =0;i<dofs_per_cell;++i)
1305 *
ad_helper.compute_linearization(cell_matrix);
1306 *
cell->get_dof_indices(local_dof_indices);
1311 *
for (
unsigned int i =0;i < dofs_per_cell; ++i){
1312 *
for(
unsigned int j = 0;
j < dofs_per_cell;++
j){
1313 *
system_matrix.add(local_dof_indices[i],local_dof_indices[
j],
cell_matrix(i,
j));
1335 *
std::cout <<
" Factorizing Jacobian matrix" << std::endl;
1343<a name=
"ann-source/compute_residual.cc"></a>
1385 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
1408 *
std::vector<types::global_dof_index> local_dof_indices(
1414 *
std::vector<double>
consol(n_q_points);
1430 *
for (
const auto &cell: dof_handler.active_cell_iterators())
1433 *
fe_values.reinit(cell);
1434 *
cell->get_dof_indices(local_dof_indices);
1435 *
const unsigned int n_independent_variables = local_dof_indices.size();
1436 *
const unsigned int n_dependent_variables = dofs_per_cell;
1445 *
fe_values[t].get_function_values_from_local_dof_values(
dof_values_ad,
1447 *
fe_values[t].get_function_gradients_from_local_dof_values(
dof_values_ad,
1458 *
std::vector<ADNumberType>
residual_ad(n_dependent_variables,
1460 *
for (
unsigned int q_index = 0; q_index < n_q_points; ++q_index)
1462 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1483 *
if (cell->face(face_number)->boundary_id() == 3)
1485 *
fe_face_values.reinit(cell, face_number);
1488 *
for (
unsigned int i =0;i<dofs_per_cell;++i)
1507 *
cell->get_dof_indices(local_dof_indices);
1509 *
for (
unsigned int i =0;i < dofs_per_cell; ++i)
1510 *
residual(local_dof_indices[i])+=
cell_rhs(i);
1517 *
std::cout <<
" The Norm is :: = " << residual.l2_norm() << std::endl;
1522<a name=
"ann-source/initial_conditions.cc"></a>
1547 *
double Initialcondition::value(
const Point<2> & ,
const unsigned int )
const
1557<a name=
"ann-source/nonlinear_heat_cons_des.cc"></a>
1598 *
dof_handler.
clear();
1605<a name=
"ann-source/output_results.cc"></a>
1631 *
data_out.attach_dof_handler(dof_handler);
1635 *
data_out.build_patches();
1639 *
data_out.write_vtu(output);
1644<a name=
"ann-source/set_boundary_conditions.cc"></a>
1684<a name=
"ann-source/setup_system.cc"></a>
1711 *
dof_handler.distribute_dofs(fe);
1718 *
sparsity_pattern.copy_from(
dsp);
1719 *
system_matrix.reinit(sparsity_pattern);
1725<a name=
"ann-source/solve_and_run.cc"></a>
1751 *
std::cout <<
" Solving linear system" << std::endl;
1755 *
void nonlinear_heat::run()
1759 *
std::ifstream f(
"mesh/mesh.msh");
1765 *
unsigned int prn =0;
1789 *
std::cout<<
">>>>> Time now is: "<<time <<std::endl;
1801 *
additional_data.max_iter = 100;
1828 *
const double tolerance) {
1829 *
solve(
rhs, dst, tolerance);
1848 *
time=time+delta_t;
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
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_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
Expression differentiate(const Expression &f, const Expression &x)
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.)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
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 > b(const Tensor< 2, dim, Number > &F)
void apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_data, const ViewTypeIn in, ViewTypeOut out)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
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)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation