|
ANS-1-1
|
00001 /* Author: Thomas Wick, 2012 */ 00002 /* University of Heidelberg */ 00003 /* Date: May 3, 2011, */ 00004 /* Revised Jan 30, 2012 */ 00005 /* E-mail: thomas.wick@iwr.uni-heidelberg.de */ 00006 /* 00007 /* */ 00008 /* Copyright (C) 2011, 2012 by Thomas Wick */ 00009 /* */ 00010 /* This file is subject to QPL and may not be distributed */ 00011 /* without copyright and license information. Please refer */ 00012 /* to the webpage XXX for the text and */ 00013 /* further information on this license. */ 00014 00015 00097 #include <base/quadrature_lib.h> 00098 #include <base/logstream.h> 00099 #include <base/function.h> 00100 #include <base/utilities.h> 00101 #include <base/timer.h> 00102 00103 #include <lac/block_vector.h> 00104 #include <lac/full_matrix.h> 00105 #include <lac/block_sparse_matrix.h> 00106 #include <lac/sparse_direct.h> 00107 #include <lac/constraint_matrix.h> 00108 00109 #include <grid/tria.h> 00110 #include <grid/grid_generator.h> 00111 #include <grid/tria_accessor.h> 00112 #include <grid/tria_iterator.h> 00113 #include <grid/tria_boundary_lib.h> 00114 #include <grid/grid_tools.h> 00115 #include <grid/grid_in.h> 00116 00117 #include <dofs/dof_handler.h> 00118 #include <dofs/dof_renumbering.h> 00119 #include <dofs/dof_accessor.h> 00120 #include <dofs/dof_tools.h> 00121 00122 #include <fe/fe_q.h> 00123 #include <fe/fe_dgq.h> 00124 #include <fe/fe_dgp.h> 00125 #include <fe/fe_system.h> 00126 #include <fe/fe_values.h> 00127 #include <fe/mapping_q1.h> 00128 00129 #include <numerics/vector_tools.h> 00130 #include <numerics/matrix_tools.h> 00131 #include <numerics/data_out.h> 00132 #include <numerics/solution_transfer.h> 00133 00134 00135 // C++ 00136 #include <fstream> 00137 #include <sstream> 00138 00139 // At the end of this top-matter, we import 00140 // all deal.II names into the global 00141 // namespace: 00142 using namespace dealii; 00143 00144 00145 00154 namespace ALETransformations 00155 { 00156 template <int dim> 00157 inline 00158 Tensor<2,dim> 00159 get_pI (unsigned int q, 00160 std::vector<Vector<double> > old_solution_values) 00161 { 00162 Tensor<2,dim> tmp; 00163 tmp[0][0] = old_solution_values[q](dim+dim); 00164 tmp[1][1] = old_solution_values[q](dim+dim); 00165 00166 return tmp; 00167 } 00168 00169 template <int dim> 00170 inline 00171 Tensor<2,dim> 00172 get_pI_LinP (const double phi_i_p) 00173 { 00174 Tensor<2,dim> tmp; 00175 tmp.clear(); 00176 tmp[0][0] = phi_i_p; 00177 tmp[1][1] = phi_i_p; 00178 00179 return tmp; 00180 } 00181 00182 template <int dim> 00183 inline 00184 Tensor<1,dim> 00185 get_grad_p (unsigned int q, 00186 std::vector<std::vector<Tensor<1,dim> > > old_solution_grads) 00187 { 00188 Tensor<1,dim> grad_p; 00189 grad_p[0] = old_solution_grads[q][dim+dim][0]; 00190 grad_p[1] = old_solution_grads[q][dim+dim][1]; 00191 00192 return grad_p; 00193 } 00194 00195 template <int dim> 00196 inline 00197 Tensor<1,dim> 00198 get_grad_p_LinP (const Tensor<1,dim> phi_i_grad_p) 00199 { 00200 Tensor<1,dim> grad_p; 00201 grad_p[0] = phi_i_grad_p[0]; 00202 grad_p[1] = phi_i_grad_p[1]; 00203 00204 return grad_p; 00205 } 00206 00207 template <int dim> 00208 inline 00209 Tensor<2,dim> 00210 get_grad_u (unsigned int q, 00211 std::vector<std::vector<Tensor<1,dim> > > old_solution_grads) 00212 { 00213 Tensor<2,dim> structure_continuation; 00214 structure_continuation[0][0] = old_solution_grads[q][dim][0]; 00215 structure_continuation[0][1] = old_solution_grads[q][dim][1]; 00216 structure_continuation[1][0] = old_solution_grads[q][dim+1][0]; 00217 structure_continuation[1][1] = old_solution_grads[q][dim+1][1]; 00218 00219 return structure_continuation; 00220 } 00221 00222 template <int dim> 00223 inline 00224 Tensor<2,dim> 00225 get_grad_v (unsigned int q, 00226 std::vector<std::vector<Tensor<1,dim> > > old_solution_grads) 00227 { 00228 Tensor<2,dim> grad_v; 00229 grad_v[0][0] = old_solution_grads[q][0][0]; 00230 grad_v[0][1] = old_solution_grads[q][0][1]; 00231 grad_v[1][0] = old_solution_grads[q][1][0]; 00232 grad_v[1][1] = old_solution_grads[q][1][1]; 00233 00234 return grad_v; 00235 } 00236 00237 template <int dim> 00238 inline 00239 Tensor<2,dim> 00240 get_grad_v_T (const Tensor<2,dim> tensor_grad_v) 00241 { 00242 Tensor<2,dim> grad_v_T; 00243 grad_v_T = transpose (tensor_grad_v); 00244 00245 return grad_v_T; 00246 } 00247 00248 template <int dim> 00249 inline 00250 Tensor<2,dim> 00251 get_grad_v_LinV (const Tensor<2,dim> phi_i_grads_v) 00252 { 00253 Tensor<2,dim> tmp; 00254 tmp[0][0] = phi_i_grads_v[0][0]; 00255 tmp[0][1] = phi_i_grads_v[0][1]; 00256 tmp[1][0] = phi_i_grads_v[1][0]; 00257 tmp[1][1] = phi_i_grads_v[1][1]; 00258 00259 return tmp; 00260 } 00261 00262 template <int dim> 00263 inline 00264 Tensor<2,dim> 00265 get_Identity () 00266 { 00267 Tensor<2,dim> identity; 00268 identity[0][0] = 1.0; 00269 identity[0][1] = 0.0; 00270 identity[1][0] = 0.0; 00271 identity[1][1] = 1.0; 00272 00273 return identity; 00274 } 00275 00276 template <int dim> 00277 inline 00278 Tensor<2,dim> 00279 get_F (unsigned int q, 00280 std::vector<std::vector<Tensor<1,dim> > > old_solution_grads) 00281 { 00282 Tensor<2,dim> F; 00283 F[0][0] = 1.0 + old_solution_grads[q][dim][0]; 00284 F[0][1] = old_solution_grads[q][dim][1]; 00285 F[1][0] = old_solution_grads[q][dim+1][0]; 00286 F[1][1] = 1.0 + old_solution_grads[q][dim+1][1]; 00287 return F; 00288 } 00289 00290 template <int dim> 00291 inline 00292 Tensor<2,dim> 00293 get_F_T (const Tensor<2,dim> F) 00294 { 00295 return transpose (F); 00296 } 00297 00298 template <int dim> 00299 inline 00300 Tensor<2,dim> 00301 get_F_Inverse (const Tensor<2,dim> F) 00302 { 00303 return invert (F); 00304 } 00305 00306 template <int dim> 00307 inline 00308 Tensor<2,dim> 00309 get_F_Inverse_T (const Tensor<2,dim> F_Inverse) 00310 { 00311 return transpose (F_Inverse); 00312 } 00313 00314 template <int dim> 00315 inline 00316 double 00317 get_J (const Tensor<2,dim> tensor_F) 00318 { 00319 return determinant (tensor_F); 00320 } 00321 00322 00323 template <int dim> 00324 inline 00325 Tensor<1,dim> 00326 get_v (unsigned int q, 00327 std::vector<Vector<double> > old_solution_values) 00328 { 00329 Tensor<1,dim> v; 00330 v[0] = old_solution_values[q](0); 00331 v[1] = old_solution_values[q](1); 00332 00333 return v; 00334 } 00335 00336 template <int dim> 00337 inline 00338 Tensor<1,dim> 00339 get_v_LinV (const Tensor<1,dim> phi_i_v) 00340 { 00341 Tensor<1,dim> tmp; 00342 tmp[0] = phi_i_v[0]; 00343 tmp[1] = phi_i_v[1]; 00344 00345 return tmp; 00346 } 00347 00348 template <int dim> 00349 inline 00350 Tensor<1,dim> 00351 get_u (unsigned int q, 00352 std::vector<Vector<double> > old_solution_values) 00353 { 00354 Tensor<1,dim> u; 00355 u[0] = old_solution_values[q](dim); 00356 u[1] = old_solution_values[q](dim+1); 00357 00358 return u; 00359 } 00360 00361 template <int dim> 00362 inline 00363 Tensor<1,dim> 00364 get_u_LinU (const Tensor<1,dim> phi_i_u) 00365 { 00366 Tensor<1,dim> tmp; 00367 tmp[0] = phi_i_u[0]; 00368 tmp[1] = phi_i_u[1]; 00369 00370 return tmp; 00371 } 00372 00373 template <int dim> 00374 inline 00375 Tensor<1,dim> 00376 get_w (unsigned int q, 00377 std::vector<Vector<double> > old_solution_values) 00378 { 00379 Tensor<1,dim> w; 00380 w[0] = old_solution_values[q](dim+dim+1); 00381 w[1] = old_solution_values[q](dim+dim+2); 00382 00383 return w; 00384 } 00385 00386 template <int dim> 00387 inline 00388 Tensor<2,dim> 00389 get_grad_w (unsigned int q, 00390 std::vector<std::vector<Tensor<1,dim> > > old_solution_grads) 00391 { 00392 Tensor<2,dim> 00393 tmp; 00394 tmp[0][0] = old_solution_grads[q][dim+dim+1][0]; 00395 tmp[0][1] = old_solution_grads[q][dim+dim+1][1]; 00396 tmp[1][0] = old_solution_grads[q][dim+dim+2][0]; 00397 tmp[1][1] = old_solution_grads[q][dim+dim+2][1]; 00398 00399 return tmp; 00400 } 00401 00402 00403 00404 template <int dim> 00405 inline 00406 double 00407 get_J_LinU (unsigned int q, 00408 const std::vector<std::vector<Tensor<1,dim> > > old_solution_grads, 00409 const Tensor<2,dim> phi_i_grads_u) 00410 { 00411 return (phi_i_grads_u[0][0] * (1 + old_solution_grads[q][dim+1][1]) + 00412 (1 + old_solution_grads[q][dim][0]) * phi_i_grads_u[1][1] - 00413 phi_i_grads_u[0][1] * old_solution_grads[q][dim+1][0] - 00414 old_solution_grads[q][dim][1] * phi_i_grads_u[1][0]); 00415 } 00416 00417 template <int dim> 00418 inline 00419 double 00420 get_J_Inverse_LinU (const double J, 00421 const double J_LinU) 00422 { 00423 return (-1.0/std::pow(J,2) * J_LinU); 00424 } 00425 00426 template <int dim> 00427 inline 00428 Tensor<2,dim> 00429 get_F_LinU (const Tensor<2,dim> phi_i_grads_u) 00430 { 00431 Tensor<2,dim> tmp; 00432 tmp[0][0] = phi_i_grads_u[0][0]; 00433 tmp[0][1] = phi_i_grads_u[0][1]; 00434 tmp[1][0] = phi_i_grads_u[1][0]; 00435 tmp[1][1] = phi_i_grads_u[1][1]; 00436 00437 return tmp; 00438 } 00439 00440 template <int dim> 00441 inline 00442 Tensor<2,dim> 00443 get_F_Inverse_LinU (const Tensor<2,dim> phi_i_grads_u, 00444 const double J, 00445 const double J_LinU, 00446 unsigned int q, 00447 std::vector<std::vector<Tensor<1,dim> > > old_solution_grads 00448 ) 00449 { 00450 Tensor<2,dim> F_tilde; 00451 F_tilde[0][0] = 1.0 + old_solution_grads[q][dim+1][1]; 00452 F_tilde[0][1] = -old_solution_grads[q][dim][1]; 00453 F_tilde[1][0] = -old_solution_grads[q][dim+1][0]; 00454 F_tilde[1][1] = 1.0 + old_solution_grads[q][dim][0]; 00455 00456 Tensor<2,dim> F_tilde_LinU; 00457 F_tilde_LinU[0][0] = phi_i_grads_u[1][1]; 00458 F_tilde_LinU[0][1] = -phi_i_grads_u[0][1]; 00459 F_tilde_LinU[1][0] = -phi_i_grads_u[1][0]; 00460 F_tilde_LinU[1][1] = phi_i_grads_u[0][0]; 00461 00462 return (-1.0/(J*J) * J_LinU * F_tilde + 00463 1.0/J * F_tilde_LinU); 00464 00465 } 00466 00467 template <int dim> 00468 inline 00469 Tensor<2,dim> 00470 get_J_F_Inverse_T_LinU (const Tensor<2,dim> phi_i_grads_u) 00471 { 00472 Tensor<2,dim> tmp; 00473 tmp[0][0] = phi_i_grads_u[1][1]; 00474 tmp[0][1] = -phi_i_grads_u[1][0]; 00475 tmp[1][0] = -phi_i_grads_u[0][1]; 00476 tmp[1][1] = phi_i_grads_u[0][0]; 00477 00478 return tmp; 00479 } 00480 00481 00482 template <int dim> 00483 inline 00484 double 00485 get_tr_C_LinU (unsigned int q, 00486 const std::vector<std::vector<Tensor<1,dim> > > old_solution_grads, 00487 const Tensor<2,dim> phi_i_grads_u) 00488 { 00489 return ((1 + old_solution_grads[q][dim][0]) * 00490 phi_i_grads_u[0][0] + 00491 old_solution_grads[q][dim][1] * 00492 phi_i_grads_u[0][1] + 00493 (1 + old_solution_grads[q][dim+1][1]) * 00494 phi_i_grads_u[1][1] + 00495 old_solution_grads[q][dim+1][0] * 00496 phi_i_grads_u[1][0]); 00497 } 00498 00499 00500 } 00501 00507 namespace NSEALE 00508 { 00509 template <int dim> 00510 inline 00511 Tensor<2,dim> 00512 get_stress_fluid_ALE (const double density, 00513 const double viscosity, 00514 const Tensor<2,dim> pI, 00515 const Tensor<2,dim> grad_v, 00516 const Tensor<2,dim> grad_v_T, 00517 const Tensor<2,dim> F_Inverse, 00518 const Tensor<2,dim> F_Inverse_T) 00519 { 00520 return (-pI + density * viscosity * 00521 (grad_v * F_Inverse + F_Inverse_T * grad_v_T )); 00522 } 00523 00524 template <int dim> 00525 inline 00526 Tensor<2,dim> 00527 get_stress_fluid_except_pressure_ALE (const double density, 00528 const double viscosity, 00529 const Tensor<2,dim> grad_v, 00530 const Tensor<2,dim> grad_v_T, 00531 const Tensor<2,dim> F_Inverse, 00532 const Tensor<2,dim> F_Inverse_T) 00533 { 00534 return (density * viscosity * (grad_v * F_Inverse + F_Inverse_T * grad_v_T)); 00535 } 00536 00537 template <int dim> 00538 inline 00539 Tensor<2,dim> 00540 get_stress_fluid_ALE_1st_term_LinAll (const Tensor<2,dim> pI, 00541 const Tensor<2,dim> F_Inverse_T, 00542 const Tensor<2,dim> J_F_Inverse_T_LinU, 00543 const Tensor<2,dim> pI_LinP, 00544 const double J) 00545 { 00546 return (-J * pI_LinP * F_Inverse_T - pI * J_F_Inverse_T_LinU); 00547 } 00548 00549 template <int dim> 00550 inline 00551 Tensor<2,dim> 00552 get_stress_fluid_ALE_2nd_term_LinAll_short (const Tensor<2,dim> J_F_Inverse_T_LinU, 00553 const Tensor<2,dim> stress_fluid_ALE, 00554 const Tensor<2,dim> grad_v, 00555 const Tensor<2,dim> grad_v_LinV, 00556 const Tensor<2,dim> F_Inverse, 00557 const Tensor<2,dim> F_Inverse_LinU, 00558 const double J, 00559 const double viscosity, 00560 const double density 00561 ) 00562 { 00563 Tensor<2,dim> sigma_LinV; 00564 Tensor<2,dim> sigma_LinU; 00565 00566 sigma_LinV = grad_v_LinV * F_Inverse + transpose(F_Inverse) * transpose(grad_v_LinV); 00567 sigma_LinU = grad_v * F_Inverse_LinU + transpose(F_Inverse_LinU) * transpose(grad_v); 00568 00569 return (density * viscosity * 00570 (sigma_LinV + sigma_LinU) * J * transpose(F_Inverse) + 00571 stress_fluid_ALE * J_F_Inverse_T_LinU); 00572 } 00573 00574 template <int dim> 00575 inline 00576 Tensor<2,dim> 00577 get_stress_fluid_ALE_3rd_term_LinAll_short (const Tensor<2,dim> F_Inverse, 00578 const Tensor<2,dim> F_Inverse_LinU, 00579 const Tensor<2,dim> grad_v, 00580 const Tensor<2,dim> grad_v_LinV, 00581 const double viscosity, 00582 const double density, 00583 const double J, 00584 const Tensor<2,dim> J_F_Inverse_T_LinU) 00585 { 00586 return density * viscosity * 00587 (J_F_Inverse_T_LinU * transpose(grad_v) * transpose(F_Inverse) + 00588 J * transpose(F_Inverse) * transpose(grad_v_LinV) * transpose(F_Inverse) + 00589 J * transpose(F_Inverse) * transpose(grad_v) * transpose(F_Inverse_LinU)); 00590 } 00591 00592 00593 00594 template <int dim> 00595 inline 00596 double 00597 get_Incompressibility_ALE (unsigned int q, 00598 std::vector<std::vector<Tensor<1,dim> > > old_solution_grads) 00599 { 00600 return (old_solution_grads[q][0][0] + 00601 old_solution_grads[q][dim+1][1] * old_solution_grads[q][0][0] - 00602 old_solution_grads[q][dim][1] * old_solution_grads[q][1][0] - 00603 old_solution_grads[q][dim+1][0] * old_solution_grads[q][0][1] + 00604 old_solution_grads[q][1][1] + 00605 old_solution_grads[q][dim][0] * old_solution_grads[q][1][1]); 00606 00607 } 00608 00609 template <int dim> 00610 inline 00611 double 00612 get_Incompressibility_ALE_LinAll (const Tensor<2,dim> phi_i_grads_v, 00613 const Tensor<2,dim> phi_i_grads_u, 00614 unsigned int q, 00615 const std::vector<std::vector<Tensor<1,dim> > > old_solution_grads) 00616 { 00617 return (phi_i_grads_v[0][0] + phi_i_grads_v[1][1] + 00618 phi_i_grads_u[1][1] * old_solution_grads[q][0][0] - 00619 phi_i_grads_u[0][1] * old_solution_grads[q][1][0] - 00620 phi_i_grads_u[1][0] * old_solution_grads[q][0][1] + 00621 phi_i_grads_u[0][0] * old_solution_grads[q][1][1]); 00622 } 00623 00624 00625 template <int dim> 00626 inline 00627 Tensor<1,dim> 00628 get_Convection_LinAll_short (const Tensor<2,dim> phi_i_grads_v, 00629 const Tensor<1,dim> phi_i_v, 00630 const double J, 00631 const double J_LinU, 00632 const Tensor<2,dim> F_Inverse, 00633 const Tensor<2,dim> F_Inverse_LinU, 00634 const Tensor<1,dim> v, 00635 const Tensor<2,dim> grad_v, 00636 const double density 00637 ) 00638 { 00639 // Linearization of fluid convection term 00640 // rho J(F^{-1}v\cdot\grad)v = rho J grad(v)F^{-1}v 00641 00642 Tensor<1,dim> convection_LinU; 00643 convection_LinU = (J_LinU * grad_v * F_Inverse * v + 00644 J * grad_v * F_Inverse_LinU * v); 00645 00646 Tensor<1,dim> convection_LinV; 00647 convection_LinV = (J * (phi_i_grads_v * F_Inverse * v + 00648 grad_v * F_Inverse * phi_i_v)); 00649 00650 return density * (convection_LinU + convection_LinV); 00651 } 00652 00653 00654 template <int dim> 00655 inline 00656 Tensor<1,dim> 00657 get_Convection_u_LinAll_short (const Tensor<2,dim> phi_i_grads_v, 00658 const Tensor<1,dim> phi_i_u, 00659 const double J, 00660 const double J_LinU, 00661 const Tensor<2,dim> F_Inverse, 00662 const Tensor<2,dim> F_Inverse_LinU, 00663 const Tensor<1,dim> u, 00664 const Tensor<2,dim> grad_v, 00665 const double density 00666 ) 00667 { 00668 // Linearization of fluid convection term 00669 // rho J(F^{-1}v\cdot\grad)u = rho J grad(v)F^{-1}u 00670 00671 Tensor<1,dim> convection_LinU; 00672 convection_LinU = (J_LinU * grad_v * F_Inverse * u + 00673 J * grad_v * F_Inverse_LinU * u + 00674 J * grad_v * F_Inverse * phi_i_u); 00675 00676 Tensor<1,dim> convection_LinV; 00677 convection_LinV = (J * phi_i_grads_v * F_Inverse * u); 00678 00679 return density * (convection_LinU + convection_LinV); 00680 } 00681 00682 00683 00684 template <int dim> 00685 inline 00686 Tensor<1,dim> 00687 get_Convection_u_old_LinAll_short (const Tensor<2,dim> phi_i_grads_v, 00688 const double J, 00689 const double J_LinU, 00690 const Tensor<2,dim> F_Inverse, 00691 const Tensor<2,dim> F_Inverse_LinU, 00692 const Tensor<1,dim> old_timestep_solution_displacement, 00693 const Tensor<2,dim> grad_v, 00694 const double density 00695 ) 00696 { 00697 // Linearization of fluid convection term 00698 // rho J(F^{-1}v\cdot\grad)u = rho J grad(v)F^{-1}u 00699 00700 Tensor<1,dim> convection_LinU; 00701 convection_LinU = (J_LinU * grad_v * F_Inverse * old_timestep_solution_displacement + 00702 J * grad_v * F_Inverse_LinU * old_timestep_solution_displacement); 00703 00704 Tensor<1,dim> convection_LinV; 00705 convection_LinV = (J * phi_i_grads_v * F_Inverse * old_timestep_solution_displacement); 00706 00707 00708 return density * (convection_LinU + convection_LinV); 00709 } 00710 00711 template <int dim> 00712 inline 00713 Tensor<1,dim> 00714 get_accelaration_term_LinAll (const Tensor<1,dim> phi_i_v, 00715 const Tensor<1,dim> v, 00716 const Tensor<1,dim> old_timestep_v, 00717 const double J_LinU, 00718 const double J, 00719 const double old_timestep_J, 00720 const double density) 00721 { 00722 return density/2.0 * (J_LinU * (v - old_timestep_v) + (J + old_timestep_J) * phi_i_v); 00723 00724 } 00725 00726 00727 } 00728 00729 00733 namespace StructureTermsALE 00734 { 00735 // Green-Lagrange strain tensor 00736 template <int dim> 00737 inline 00738 Tensor<2,dim> 00739 get_E (const Tensor<2,dim> F_T, 00740 const Tensor<2,dim> F, 00741 const Tensor<2,dim> Identity) 00742 { 00743 return 0.5 * (F_T * F - Identity); 00744 } 00745 00746 template <int dim> 00747 inline 00748 double 00749 get_tr_E (const Tensor<2,dim> E) 00750 { 00751 return trace (E); 00752 } 00753 00754 template <int dim> 00755 inline 00756 double 00757 get_tr_E_LinU (unsigned int q, 00758 const std::vector<std::vector<Tensor<1,dim> > > old_solution_grads, 00759 const Tensor<2,dim> phi_i_grads_u) 00760 { 00761 return ((1 + old_solution_grads[q][dim][0]) * 00762 phi_i_grads_u[0][0] + 00763 old_solution_grads[q][dim][1] * 00764 phi_i_grads_u[0][1] + 00765 (1 + old_solution_grads[q][dim+1][1]) * 00766 phi_i_grads_u[1][1] + 00767 old_solution_grads[q][dim+1][0] * 00768 phi_i_grads_u[1][0]); 00769 } 00770 00771 } 00772 00773 00774 00784 template <int dim> 00785 class BoundaryParabolic : public Function<dim> 00786 { 00787 public: 00788 BoundaryParabolic (const double time) 00789 : Function<dim>(dim+dim+dim+1) 00790 { 00791 _time = time; 00792 } 00793 00794 virtual double value (const Point<dim> &p, 00795 const unsigned int component = 0) const; 00796 00797 virtual void vector_value (const Point<dim> &p, 00798 Vector<double> &value) const; 00799 00800 private: 00801 double _time; 00802 00803 }; 00804 00805 // The boundary values are given to component 00806 // with number 0. 00807 template <int dim> 00808 double 00809 BoundaryParabolic<dim>::value (const Point<dim> &p, 00810 const unsigned int component) const 00811 { 00812 Assert (component < this->n_components, 00813 ExcIndexRange (component, 0, this->n_components)); 00814 00815 const long double pi = 3.141592653589793238462643; 00816 00817 // The maximum inflow depends on the configuration 00818 // for the different test cases: 00819 // FSI 1: 0.2; FSI 2: 1.0; FSI 3: 2.0 00820 // 00821 // For the two unsteady test cases FSI 2 and FSI 3, it 00822 // is recommanded to start with a smooth increase of 00823 // the inflow. Hence, we use the cosine function 00824 // to control the inflow at the beginning until 00825 // the total time 2.0 has been reached. 00826 double inflow_velocity = 2.0e-1; 00827 00828 if (component == 0) 00829 { 00830 if (_time < 2.0) 00831 { 00832 return ( (p(0) == 0) && (p(1) <= 0.41) ? -1.5 * inflow_velocity * 00833 (1.0 - std::cos(pi/2.0 * _time))/2.0 * 00834 (4.0/0.1681) * 00835 (std::pow(p(1), 2) - 0.41 * std::pow(p(1),1)) : 0 ); 00836 } 00837 else 00838 { 00839 return ( (p(0) == 0) && (p(1) <= 0.41) ? -1.5 * inflow_velocity * 00840 (4.0/0.1681) * 00841 (std::pow(p(1), 2) - 0.41 * std::pow(p(1),1)) : 0 ); 00842 00843 } 00844 00845 } 00846 00847 return 0; 00848 } 00849 00850 00851 00852 template <int dim> 00853 void 00854 BoundaryParabolic<dim>::vector_value (const Point<dim> &p, 00855 Vector<double> &values) const 00856 { 00857 for (unsigned int c=0; c<this->n_components; ++c) 00858 values (c) = BoundaryParabolic<dim>::value (p, c); 00859 } 00860 00861 00913 template <int dim> 00914 class FSIALEProblem 00915 { 00916 public: 00917 00918 FSIALEProblem (const unsigned int degree); 00919 ~FSIALEProblem (); 00920 void run (); 00921 00922 private: 00923 00924 void set_runtime_parameters (); 00925 void setup_system (); 00926 void assemble_system_matrix (); 00927 void assemble_system_rhs (); 00928 00929 void set_initial_bc (const double time); 00930 void set_newton_bc (); 00931 00932 void solve (); 00933 void newton_iteration(const double time); 00934 void output_results (const unsigned int refinement_cycle, 00935 const BlockVector<double> solution) const; 00936 00937 double compute_point_value (Point<dim> p, 00938 const unsigned int component) const; 00939 00940 void compute_drag_lift_fsi_fluid_tensor (); 00941 void compute_functional_values (); 00942 00943 const unsigned int degree; 00944 00945 Triangulation<dim> triangulation; 00946 FESystem<dim> fe; 00947 DoFHandler<dim> dof_handler; 00948 00949 ConstraintMatrix constraints; 00950 00951 BlockSparsityPattern sparsity_pattern; 00952 BlockSparseMatrix<double> system_matrix; 00953 00954 BlockVector<double> solution, newton_update, old_timestep_solution; 00955 BlockVector<double> system_rhs; 00956 00957 TimerOutput timer; 00958 00959 // Global variables for timestepping scheme 00960 unsigned int timestep_number; 00961 unsigned int max_no_timesteps; 00962 double timestep, theta, time; 00963 std::string time_stepping_scheme; 00964 00965 // Fluid parameters 00966 double density_fluid, viscosity; 00967 00968 // Structure parameters 00969 double density_structure; 00970 double lame_coefficient_mu, lame_coefficient_lambda, poisson_ratio_nu; 00971 00972 // Other parameters to control the fluid mesh motion 00973 double cell_diameter; 00974 double alpha_u, alpha_w; 00975 00976 00977 00978 00979 00980 }; 00981 00982 00989 template <int dim> 00990 FSIALEProblem<dim>::FSIALEProblem (const unsigned int degree) 00991 : 00992 degree (degree), 00993 triangulation (Triangulation<dim>::maximum_smoothing), 00994 fe (FE_Q<dim>(degree+1), dim, 00995 FE_Q<dim>(degree+1), dim, 00996 FE_DGP<dim>(degree), 1, 00997 FE_Q<dim>(degree+1), dim), 00998 dof_handler (triangulation), 00999 timer (std::cout, TimerOutput::summary, TimerOutput::cpu_times) 01000 {} 01001 01002 01003 // This is the standard destructor. 01004 template <int dim> 01005 FSIALEProblem<dim>::~FSIALEProblem () 01006 {} 01007 01008 01015 template <int dim> 01016 void FSIALEProblem<dim>::set_runtime_parameters () 01017 { 01018 // Fluid parameters 01019 density_fluid = 1.0e+3; 01020 01021 // FSI 1 & 3: 1.0e+3; FSI 2: 1.0e+4 01022 density_structure = 1.0e+3; 01023 viscosity = 1.0e-3; 01024 01025 // Structure parameters 01026 // FSI 1 & 2: 0.5e+6; FSI 3: 2.0e+6 01027 lame_coefficient_mu = 0.5e+6; 01028 poisson_ratio_nu = 0.4; 01029 01030 lame_coefficient_lambda = (2 * poisson_ratio_nu * lame_coefficient_mu)/ 01031 (1.0 - 2 * poisson_ratio_nu); 01032 01033 // Diffusion parameters to control the fluid mesh motion 01034 // The higher these parameters the stiffer the fluid mesh. 01035 alpha_u = 1.0e-5; 01036 alpha_w = 1.0e-5; 01037 01038 // Timestepping schemes 01039 //BE, CN, CN_shifted 01040 time_stepping_scheme = "BE"; 01041 01042 // Timestep size: 01043 // FSI 1: 1.0 (quasi-stationary) 01044 // FSI 2: <= 1.0e-2 (non-stationary) 01045 // FSI 3: <= 1.0e-3 (non-stationary) 01046 timestep = 1.0; 01047 01048 // Maximum number of timesteps: 01049 // FSI 1: 25 , T= 25 (timestep == 1.0) 01050 // FSI 2: 1500, T= 15 (timestep == 1.0e-2) 01051 // FSI 3: 10000, T= 10 (timestep == 1.0e-3) 01052 max_no_timesteps = 25; 01053 01054 // A variable to count the number of time steps 01055 timestep_number = 0; 01056 01057 // Counts total time 01058 time = 0; 01059 01060 // Here, we choose a time-stepping scheme that 01061 // is based on finite differences: 01062 // BE = backward Euler scheme 01063 // CN = Crank-Nicolson scheme 01064 // CN_shifted = time-shifted Crank-Nicolson scheme 01065 // For further properties of these schemes, 01066 // we refer to standard literature. 01067 if (time_stepping_scheme == "BE") 01068 theta = 1.0; 01069 else if (time_stepping_scheme == "CN") 01070 theta = 0.5; 01071 else if (time_stepping_scheme == "CN_shifted") 01072 theta = 0.5 + timestep; 01073 else 01074 std::cout << "No such timestepping scheme" << std::endl; 01075 01076 // In the following, we read a *.inp grid from a file. 01077 // The geometry information is based on the 01078 // fluid-structure interaction benchmark problems 01079 // (Lit. J. Hron, S. Turek, 2006) 01080 std::string grid_name; 01081 grid_name = "fsi.inp"; 01082 01083 GridIn<dim> grid_in; 01084 grid_in.attach_triangulation (triangulation); 01085 std::ifstream input_file(grid_name.c_str()); 01086 Assert (dim==2, ExcInternalError()); 01087 grid_in.read_ucd (input_file); 01088 01089 Point<dim> p(0.2, 0.2); 01090 double radius = 0.05; 01091 static const HyperBallBoundary<dim> boundary(p,radius); 01092 triangulation.set_boundary (80, boundary); 01093 triangulation.set_boundary (81, boundary); 01094 01095 triangulation.refine_global (1); 01096 01097 } 01098 01099 01100 01101 // This function is similar to many deal.II tuturial steps. 01102 template <int dim> 01103 void FSIALEProblem<dim>::setup_system () 01104 { 01105 timer.enter_section("Setup system."); 01106 01107 // We set runtime parameters to drive the problem. 01108 // These parameters could also be read from a parameter file that 01109 // can be handled by the ParameterHandler object (see step-19) 01110 set_runtime_parameters (); 01111 01112 system_matrix.clear (); 01113 01114 dof_handler.distribute_dofs (fe); 01115 DoFRenumbering::Cuthill_McKee (dof_handler); 01116 01117 // We are dealing with 7 components for this 01118 // two-dimensional fluid-structure interacion problem 01119 // Precisely, we use: 01120 // velocity in x and y: 0 01121 // structure displacement in x and y: 1 01122 // scalar pressure field: 2 01123 // additional displacement in x and y: 3 01124 std::vector<unsigned int> block_component (7,0); 01125 block_component[dim] = 1; 01126 block_component[dim+1] = 1; 01127 block_component[dim+dim] = 2; 01128 block_component[dim+dim+1] = 3; 01129 block_component[dim+dim+dim] = 3; 01130 01131 DoFRenumbering::component_wise (dof_handler, block_component); 01132 01133 { 01134 constraints.clear (); 01135 set_newton_bc (); 01136 DoFTools::make_hanging_node_constraints (dof_handler, 01137 constraints); 01138 } 01139 constraints.close (); 01140 01141 std::vector<unsigned int> dofs_per_block (4); 01142 DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component); 01143 const unsigned int n_v = dofs_per_block[0], 01144 n_u = dofs_per_block[1], 01145 n_p = dofs_per_block[2], 01146 n_w = dofs_per_block[3]; 01147 01148 std::cout << "Cells:\t" 01149 << triangulation.n_active_cells() 01150 << std::endl 01151 << "DoFs:\t" 01152 << dof_handler.n_dofs() 01153 << " (" << n_v << '+' << n_u << '+' << n_p << '+' << n_w << ')' 01154 << std::endl; 01155 01156 01157 01158 01159 { 01160 BlockCompressedSimpleSparsityPattern csp (4,4); 01161 01162 csp.block(0,0).reinit (n_v, n_v); 01163 csp.block(0,1).reinit (n_v, n_u); 01164 csp.block(0,2).reinit (n_v, n_p); 01165 csp.block(0,3).reinit (n_v, n_w); 01166 01167 csp.block(1,0).reinit (n_u, n_v); 01168 csp.block(1,1).reinit (n_u, n_u); 01169 csp.block(1,2).reinit (n_u, n_p); 01170 csp.block(1,3).reinit (n_u, n_w); 01171 01172 csp.block(2,0).reinit (n_p, n_v); 01173 csp.block(2,1).reinit (n_p, n_u); 01174 csp.block(2,2).reinit (n_p, n_p); 01175 csp.block(2,3).reinit (n_p, n_w); 01176 01177 csp.block(3,0).reinit (n_w, n_v); 01178 csp.block(3,1).reinit (n_w, n_u); 01179 csp.block(3,2).reinit (n_w, n_p); 01180 csp.block(3,3).reinit (n_w, n_w); 01181 01182 csp.collect_sizes(); 01183 01184 01185 DoFTools::make_sparsity_pattern (dof_handler, csp, constraints, false); 01186 01187 sparsity_pattern.copy_from (csp); 01188 } 01189 01190 system_matrix.reinit (sparsity_pattern); 01191 01192 // Actual solution at time step n 01193 solution.reinit (4); 01194 solution.block(0).reinit (n_v); 01195 solution.block(1).reinit (n_u); 01196 solution.block(2).reinit (n_p); 01197 solution.block(3).reinit (n_w); 01198 01199 solution.collect_sizes (); 01200 01201 // Old timestep solution at time step n-1 01202 old_timestep_solution.reinit (4); 01203 old_timestep_solution.block(0).reinit (n_v); 01204 old_timestep_solution.block(1).reinit (n_u); 01205 old_timestep_solution.block(2).reinit (n_p); 01206 old_timestep_solution.block(3).reinit (n_w); 01207 01208 old_timestep_solution.collect_sizes (); 01209 01210 01211 // Updates for Newton's method 01212 newton_update.reinit (4); 01213 newton_update.block(0).reinit (n_v); 01214 newton_update.block(1).reinit (n_u); 01215 newton_update.block(2).reinit (n_p); 01216 newton_update.block(3).reinit (n_w); 01217 01218 newton_update.collect_sizes (); 01219 01220 // Residual for Newton's method 01221 system_rhs.reinit (4); 01222 system_rhs.block(0).reinit (n_v); 01223 system_rhs.block(1).reinit (n_u); 01224 system_rhs.block(2).reinit (n_p); 01225 system_rhs.block(3).reinit (n_w); 01226 01227 system_rhs.collect_sizes (); 01228 01229 timer.exit_section(); 01230 } 01231 01232 01258 template <int dim> 01259 void FSIALEProblem<dim>::assemble_system_matrix () 01260 { 01261 timer.enter_section("Assemble Matrix."); 01262 system_matrix=0; 01263 01264 QGauss<dim> quadrature_formula(degree+2); 01265 QGauss<dim-1> face_quadrature_formula(degree+2); 01266 01267 FEValues<dim> fe_values (fe, quadrature_formula, 01268 update_values | 01269 update_quadrature_points | 01270 update_JxW_values | 01271 update_gradients); 01272 01273 FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula, 01274 update_values | update_quadrature_points | 01275 update_normal_vectors | update_gradients | 01276 update_JxW_values); 01277 01278 const unsigned int dofs_per_cell = fe.dofs_per_cell; 01279 01280 const unsigned int n_q_points = quadrature_formula.size(); 01281 const unsigned int n_face_q_points = face_quadrature_formula.size(); 01282 01283 FullMatrix<double> local_matrix (dofs_per_cell, dofs_per_cell); 01284 01285 std::vector<unsigned int> local_dof_indices (dofs_per_cell); 01286 01287 01288 // Now, we are going to use the 01289 // FEValuesExtractors to determine 01290 // the four principle variables 01291 const FEValuesExtractors::Vector velocities (0); 01292 const FEValuesExtractors::Vector displacements (dim); // 2 01293 const FEValuesExtractors::Scalar pressure (dim+dim); // 4 01294 const FEValuesExtractors::Vector displacements_w (dim+dim+1); // 4 01295 01296 01297 // We declare Vectors and Tensors for 01298 // the solutions at the previous Newton iteration: 01299 std::vector<Vector<double> > old_solution_values (n_q_points, 01300 Vector<double>(dim+dim+dim+1)); 01301 01302 std::vector<std::vector<Tensor<1,dim> > > old_solution_grads (n_q_points, 01303 std::vector<Tensor<1,dim> > (dim+dim+dim+1)); 01304 01305 std::vector<Vector<double> > old_solution_face_values (n_face_q_points, 01306 Vector<double>(dim+dim+dim+1)); 01307 01308 std::vector<std::vector<Tensor<1,dim> > > old_solution_face_grads (n_face_q_points, 01309 std::vector<Tensor<1,dim> > (dim+dim+dim+1)); 01310 01311 01312 // We declare Vectors and Tensors for 01313 // the solution at the previous time step: 01314 std::vector<Vector<double> > old_timestep_solution_values (n_q_points, 01315 Vector<double>(dim+dim+dim+1)); 01316 01317 01318 std::vector<std::vector<Tensor<1,dim> > > old_timestep_solution_grads (n_q_points, 01319 std::vector<Tensor<1,dim> > (dim+dim+dim+1)); 01320 01321 01322 std::vector<Vector<double> > old_timestep_solution_face_values (n_face_q_points, 01323 Vector<double>(dim+dim+dim+1)); 01324 01325 01326 std::vector<std::vector<Tensor<1,dim> > > old_timestep_solution_face_grads (n_face_q_points, 01327 std::vector<Tensor<1,dim> > (dim+dim+dim+1)); 01328 01329 // Declaring test functions: 01330 std::vector<Tensor<1,dim> > phi_i_v (dofs_per_cell); 01331 std::vector<Tensor<2,dim> > phi_i_grads_v(dofs_per_cell); 01332 std::vector<double> phi_i_p(dofs_per_cell); 01333 std::vector<Tensor<1,dim> > phi_i_u (dofs_per_cell); 01334 std::vector<Tensor<2,dim> > phi_i_grads_u(dofs_per_cell); 01335 std::vector<Tensor<1,dim> > phi_i_w (dofs_per_cell); 01336 std::vector<Tensor<2,dim> > phi_i_grads_w(dofs_per_cell); 01337 01338 // This is the identity matrix in two dimensions: 01339 const Tensor<2,dim> Identity = ALETransformations 01340 ::get_Identity<dim> (); 01341 01342 typename DoFHandler<dim>::active_cell_iterator 01343 cell = dof_handler.begin_active(), 01344 endc = dof_handler.end(); 01345 01346 for (; cell!=endc; ++cell) 01347 { 01348 fe_values.reinit (cell); 01349 local_matrix = 0; 01350 01351 // We need the cell diameter to control the fluid mesh motion 01352 cell_diameter = cell->diameter(); 01353 01354 // Old Newton iteration values 01355 fe_values.get_function_values (solution, old_solution_values); 01356 fe_values.get_function_grads (solution, old_solution_grads); 01357 01358 // Old_timestep_solution values 01359 fe_values.get_function_values (old_timestep_solution, old_timestep_solution_values); 01360 fe_values.get_function_grads (old_timestep_solution, old_timestep_solution_grads); 01361 01362 // Next, we run over all cells for the fluid equations 01363 if (cell->material_id() == 0) 01364 { 01365 for (unsigned int q=0; q<n_q_points; ++q) 01366 { 01367 for (unsigned int k=0; k<dofs_per_cell; ++k) 01368 { 01369 phi_i_v[k] = fe_values[velocities].value (k, q); 01370 phi_i_grads_v[k] = fe_values[velocities].gradient (k, q); 01371 phi_i_p[k] = fe_values[pressure].value (k, q); 01372 phi_i_u[k] = fe_values[displacements].value (k, q); 01373 phi_i_grads_u[k] = fe_values[displacements].gradient (k, q); 01374 phi_i_w[k] = fe_values[displacements_w].value (k, q); 01375 phi_i_grads_w[k] = fe_values[displacements_w].gradient (k, q); 01376 } 01377 01378 // We build values, vectors, and tensors 01379 // from information of the previous Newton step. These are introduced 01380 // for two reasons: 01381 // First, these are used to perform the ALE mapping of the 01382 // fluid equations. Second, these terms are used to 01383 // make the notation as simple and self-explaining as possible: 01384 const Tensor<2,dim> pI = ALETransformations 01385 ::get_pI<dim> (q, old_solution_values); 01386 01387 const Tensor<1,dim> v = ALETransformations 01388 ::get_v<dim> (q, old_solution_values); 01389 01390 const Tensor<1,dim> u = ALETransformations 01391 ::get_u<dim> (q,old_solution_values); 01392 01393 const Tensor<2,dim> grad_v = ALETransformations 01394 ::get_grad_v<dim> (q, old_solution_grads); 01395 01396 const Tensor<2,dim> grad_v_T = ALETransformations 01397 ::get_grad_v_T<dim> (grad_v); 01398 01399 const Tensor<2,dim> F = ALETransformations 01400 ::get_F<dim> (q, old_solution_grads); 01401 01402 const Tensor<2,dim> F_Inverse = ALETransformations 01403 ::get_F_Inverse<dim> (F); 01404 01405 const Tensor<2,dim> F_Inverse_T = ALETransformations 01406 ::get_F_Inverse_T<dim> (F_Inverse); 01407 01408 const double J = ALETransformations 01409 ::get_J<dim> (F); 01410 01411 // Stress tensor for the fluid in ALE notation 01412 const Tensor<2,dim> sigma_ALE = NSEALE 01413 ::get_stress_fluid_ALE<dim> (density_fluid, viscosity, pI, 01414 grad_v, grad_v_T, F_Inverse, F_Inverse_T ); 01415 01416 // Further, we also need some information from the previous time steps 01417 const Tensor<1,dim> old_timestep_v = ALETransformations 01418 ::get_v<dim> (q, old_timestep_solution_values); 01419 01420 const Tensor<1,dim> old_timestep_u = ALETransformations 01421 ::get_u<dim> (q, old_timestep_solution_values); 01422 01423 const Tensor<2,dim> old_timestep_F = ALETransformations 01424 ::get_F<dim> (q, old_timestep_solution_grads); 01425 01426 const double old_timestep_J = ALETransformations 01427 ::get_J<dim> (old_timestep_F); 01428 01429 // Outer loop for dofs 01430 for (unsigned int i=0; i<dofs_per_cell; ++i) 01431 { 01432 const Tensor<2,dim> pI_LinP = ALETransformations 01433 ::get_pI_LinP<dim> (phi_i_p[i]); 01434 01435 const Tensor<2,dim> grad_v_LinV = ALETransformations 01436 ::get_grad_v_LinV<dim> (phi_i_grads_v[i]); 01437 01438 const double J_LinU = ALETransformations 01439 ::get_J_LinU<dim> (q, old_solution_grads, phi_i_grads_u[i]); 01440 01441 const Tensor<2,dim> J_F_Inverse_T_LinU = ALETransformations 01442 ::get_J_F_Inverse_T_LinU<dim> (phi_i_grads_u[i]); 01443 01444 const Tensor<2,dim> F_Inverse_LinU = ALETransformations 01445 ::get_F_Inverse_LinU (phi_i_grads_u[i], J, J_LinU, q, old_solution_grads); 01446 01447 const Tensor<2,dim> stress_fluid_ALE_1st_term_LinAll = NSEALE 01448 ::get_stress_fluid_ALE_1st_term_LinAll<dim> 01449 (pI, F_Inverse_T, J_F_Inverse_T_LinU, pI_LinP, J); 01450 01451 const Tensor<2,dim> stress_fluid_ALE_2nd_term_LinAll = NSEALE 01452 ::get_stress_fluid_ALE_2nd_term_LinAll_short 01453 (J_F_Inverse_T_LinU, sigma_ALE, grad_v, grad_v_LinV, 01454 F_Inverse, F_Inverse_LinU, J, viscosity, density_fluid); 01455 01456 const Tensor<1,dim> convection_fluid_LinAll_short = NSEALE 01457 ::get_Convection_LinAll_short<dim> 01458 (phi_i_grads_v[i], phi_i_v[i], J,J_LinU, 01459 F_Inverse, F_Inverse_LinU, v, grad_v, density_fluid); 01460 01461 const double incompressibility_ALE_LinAll = NSEALE 01462 ::get_Incompressibility_ALE_LinAll<dim> 01463 (phi_i_grads_v[i], phi_i_grads_u[i], q, old_solution_grads); 01464 01465 const Tensor<1,dim> accelaration_term_LinAll = NSEALE 01466 ::get_accelaration_term_LinAll 01467 (phi_i_v[i], v, old_timestep_v, J_LinU, 01468 J, old_timestep_J, density_fluid); 01469 01470 const Tensor<1,dim> convection_fluid_u_LinAll_short = NSEALE 01471 ::get_Convection_u_LinAll_short<dim> 01472 (phi_i_grads_v[i], phi_i_u[i], J,J_LinU, F_Inverse, 01473 F_Inverse_LinU, u, grad_v, density_fluid); 01474 01475 const Tensor<1,dim> convection_fluid_u_old_LinAll_short = NSEALE 01476 ::get_Convection_u_old_LinAll_short<dim> 01477 (phi_i_grads_v[i], J, J_LinU, F_Inverse, 01478 F_Inverse_LinU, old_timestep_u, grad_v, density_fluid); 01479 01480 // Inner loop for dofs 01481 for (unsigned int j=0; j<dofs_per_cell; ++j) 01482 { 01483 // Fluid , NSE in ALE 01484 const unsigned int comp_j = fe.system_to_component_index(j).first; 01485 if (comp_j == 0 || comp_j == 1) 01486 { 01487 local_matrix(j,i) += (accelaration_term_LinAll * phi_i_v[j] + 01488 timestep * theta * 01489 convection_fluid_LinAll_short * phi_i_v[j] - 01490 convection_fluid_u_LinAll_short * phi_i_v[j] + 01491 convection_fluid_u_old_LinAll_short * phi_i_v[j] + 01492 timestep * scalar_product(stress_fluid_ALE_1st_term_LinAll, phi_i_grads_v[j]) + 01493 timestep * theta * 01494 scalar_product(stress_fluid_ALE_2nd_term_LinAll, phi_i_grads_v[j]) 01495 ) * fe_values.JxW(q); 01496 } 01497 else if (comp_j == 2 || comp_j == 3) 01498 { 01499 local_matrix(j,i) += (alpha_u * scalar_product(phi_i_grads_w[i], phi_i_grads_u[j]) 01500 ) * fe_values.JxW(q); 01501 } 01502 else if (comp_j == 4) 01503 { 01504 local_matrix(j,i) += (incompressibility_ALE_LinAll * phi_i_p[j] 01505 ) * fe_values.JxW(q); 01506 } 01507 else if (comp_j == 5 || comp_j == 6) 01508 { 01509 local_matrix(j,i) += (alpha_w * (phi_i_w[i] * phi_i_w[j] - scalar_product(phi_i_grads_u[i],phi_i_grads_w[j])) 01510 ) * fe_values.JxW(q); 01511 } 01512 01513 // end j dofs 01514 } 01515 // end i dofs 01516 } 01517 // end n_q_points 01518 } 01519 01520 // We compute in the following 01521 // one term on the outflow boundary. 01522 // This relation is well-know in the literature 01523 // as "do-nothing" condition. Therefore, we only 01524 // ask for the corresponding color at the outflow 01525 // boundary that is 1 in our case. 01526 for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face) 01527 { 01528 if (cell->face(face)->at_boundary() && 01529 (cell->face(face)->boundary_indicator() == 1) 01530 ) 01531 { 01532 01533 fe_face_values.reinit (cell, face); 01534 01535 fe_face_values.get_function_values (solution, old_solution_face_values); 01536 fe_face_values.get_function_grads (solution, old_solution_face_grads); 01537 01538 for (unsigned int q=0; q<n_face_q_points; ++q) 01539 { 01540 for (unsigned int k=0; k<dofs_per_cell; ++k) 01541 { 01542 phi_i_v[k] = fe_face_values[velocities].value (k, q); 01543 phi_i_grads_v[k] = fe_face_values[velocities].gradient (k, q); 01544 phi_i_grads_u[k] = fe_face_values[displacements].gradient (k, q); 01545 } 01546 01547 const Tensor<2,dim> pI = ALETransformations 01548 ::get_pI<dim> (q, old_solution_face_values); 01549 01550 const Tensor<1,dim> v = ALETransformations 01551 ::get_v<dim> (q, old_solution_face_values); 01552 01553 const Tensor<2,dim> grad_v = ALETransformations 01554 ::get_grad_v<dim> (q, old_solution_face_grads); 01555 01556 const Tensor<2,dim> grad_v_T = ALETransformations 01557 ::get_grad_v_T<dim> (grad_v); 01558 01559 const Tensor<2,dim> F = ALETransformations 01560 ::get_F<dim> (q, old_solution_face_grads); 01561 01562 const Tensor<2,dim> F_Inverse = ALETransformations 01563 ::get_F_Inverse<dim> (F); 01564 01565 const Tensor<2,dim> F_Inverse_T = ALETransformations 01566 ::get_F_Inverse_T<dim> (F_Inverse); 01567 01568 const double J = ALETransformations 01569 ::get_J<dim> (F); 01570 01571 01572 for (unsigned int i=0; i<dofs_per_cell; ++i) 01573 { 01574 const Tensor<2,dim> grad_v_LinV = ALETransformations 01575 ::get_grad_v_LinV<dim> (phi_i_grads_v[i]); 01576 01577 const double J_LinU = ALETransformations 01578 ::get_J_LinU<dim> (q, old_solution_face_grads, phi_i_grads_u[i]); 01579 01580 const Tensor<2,dim> J_F_Inverse_T_LinU = ALETransformations 01581 ::get_J_F_Inverse_T_LinU<dim> (phi_i_grads_u[i]); 01582 01583 const Tensor<2,dim> F_Inverse_LinU = ALETransformations 01584 ::get_F_Inverse_LinU 01585 (phi_i_grads_u[i], J, J_LinU, q, old_solution_face_grads); 01586 01587 const Tensor<2,dim> stress_fluid_ALE_3rd_term_LinAll = NSEALE 01588 ::get_stress_fluid_ALE_3rd_term_LinAll_short<dim> 01589 (F_Inverse, F_Inverse_LinU, grad_v, grad_v_LinV, 01590 viscosity, density_fluid, J, J_F_Inverse_T_LinU); 01591 01592 // Here, we multiply the symmetric part of fluid's stress tensor 01593 // with the normal direction. 01594 const Tensor<1,dim> neumann_value 01595 = (stress_fluid_ALE_3rd_term_LinAll * fe_face_values.normal_vector(q)); 01596 01597 for (unsigned int j=0; j<dofs_per_cell; ++j) 01598 { 01599 const unsigned int comp_j = fe.system_to_component_index(j).first; 01600 if (comp_j == 0 || comp_j == 1) 01601 { 01602 local_matrix(j,i) -= (timestep * theta * 01603 neumann_value * phi_i_v[j] 01604 ) * fe_face_values.JxW(q); 01605 } 01606 // end j 01607 } 01608 // end i 01609 } 01610 // end q_face_points 01611 } 01612 // end if-routine face integrals 01613 } 01614 // end face integrals 01615 } 01616 01617 01618 // Next, we compute the face integrals on the interface between fluid and structure. 01619 // The appear because of partial integration of the additional equations for 01620 // the fluid mesh motion. 01621 for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face) 01622 { 01623 if (cell->neighbor_index(face) != -1) 01624 if (cell->material_id() != cell->neighbor(face)->material_id() && 01625 cell->face(face)->boundary_indicator()!=81) 01626 { 01627 01628 fe_face_values.reinit (cell, face); 01629 01630 for (unsigned int q=0; q<n_face_q_points; ++q) 01631 { 01632 for (unsigned int k=0; k<dofs_per_cell; ++k) 01633 { 01634 phi_i_u[k] = fe_face_values[displacements].value (k, q); 01635 phi_i_grads_u[k] = fe_face_values[displacements].gradient (k, q); 01636 01637 phi_i_w[k] = fe_face_values[displacements_w].value (k, q); 01638 phi_i_grads_w[k] = fe_face_values[displacements_w].gradient (k, q); 01639 } 01640 01641 for (unsigned int i=0; i<dofs_per_cell; ++i) 01642 { 01643 const Tensor<1,dim> neumann_value_w 01644 = (phi_i_grads_w[i] * fe_face_values.normal_vector(q)); 01645 01646 const Tensor<1,dim> neumann_value_u 01647 = (phi_i_grads_u[i] * fe_face_values.normal_vector(q)); 01648 01649 for (unsigned int j=0; j<dofs_per_cell; ++j) 01650 { 01651 const unsigned int comp_j = fe.system_to_component_index(j).first; 01652 if (comp_j == 2 || comp_j == 3) 01653 { 01654 local_matrix(j,i) -= (alpha_u * neumann_value_w * phi_i_u[j] * fe_face_values.JxW(q)); 01655 01656 } 01657 else if (comp_j == 5 || comp_j == 6) 01658 { 01659 local_matrix(j,i) -= (alpha_w * neumann_value_u * phi_i_w[j] * fe_face_values.JxW(q)); 01660 01661 } 01662 01663 } // end j 01664 01665 } // end i 01666 01667 } // end q_face_points 01668 } // end if-routine face integrals 01669 01670 01671 } // end face artificial Gamma i // 01672 01673 01674 // This is the same as discussed in step-22: 01675 cell->get_dof_indices (local_dof_indices); 01676 constraints.distribute_local_to_global (local_matrix, local_dof_indices, 01677 system_matrix); 01678 01679 // Finally, we arrive at the end for assembling the matrix 01680 // for the fluid equations and step to the computation of the 01681 // structure terms: 01682 } 01683 else if (cell->material_id() == 1) 01684 { 01685 for (unsigned int q=0; q<n_q_points; ++q) 01686 { 01687 for (unsigned int k=0; k<dofs_per_cell; ++k) 01688 { 01689 phi_i_v[k] = fe_values[velocities].value (k, q); 01690 phi_i_grads_v[k] = fe_values[velocities].gradient (k, q); 01691 phi_i_p[k] = fe_values[pressure].value (k, q); 01692 phi_i_u[k] = fe_values[displacements].value (k, q); 01693 phi_i_grads_u[k] = fe_values[displacements].gradient (k, q); 01694 phi_i_w[k] = fe_values[displacements_w].value (k, q); 01695 phi_i_grads_w[k] = fe_values[displacements_w].gradient (k, q); 01696 } 01697 01698 // It is here the same as already shown for the fluid equations. 01699 // First, we prepare things coming from the previous Newton 01700 // iteration... 01701 const Tensor<2,dim> pI = ALETransformations 01702 ::get_pI<dim> (q, old_solution_values); 01703 01704 const Tensor<2,dim> grad_v = ALETransformations 01705 ::get_grad_v<dim> (q, old_solution_grads); 01706 01707 const Tensor<2,dim> grad_v_T = ALETransformations 01708 ::get_grad_v_T<dim> (grad_v); 01709 01710 const Tensor<1,dim> v = ALETransformations 01711 ::get_v<dim> (q, old_solution_values); 01712 01713 const Tensor<2,dim> F = ALETransformations 01714 ::get_F<dim> (q, old_solution_grads); 01715 01716 const Tensor<2,dim> F_Inverse = ALETransformations 01717 ::get_F_Inverse<dim> (F); 01718 01719 const Tensor<2,dim> F_Inverse_T = ALETransformations 01720 ::get_F_Inverse_T<dim> (F_Inverse); 01721 01722 const Tensor<2,dim> F_T = ALETransformations 01723 ::get_F_T<dim> (F); 01724 01725 const double J = ALETransformations 01726 ::get_J<dim> (F); 01727 01728 const Tensor<2,dim> E = StructureTermsALE 01729 ::get_E<dim> (F_T, F, Identity); 01730 01731 const double tr_E = StructureTermsALE 01732 ::get_tr_E<dim> (E); 01733 01734 // ... and then things coming from the previous time steps 01735 const Tensor<1,dim> old_timestep_v = ALETransformations 01736 ::get_v<dim> (q, old_timestep_solution_values); 01737 01738 const Tensor<2,dim> old_timestep_F = ALETransformations 01739 ::get_F<dim> (q, old_timestep_solution_grads); 01740 01741 const double old_timestep_J = ALETransformations 01742 ::get_J<dim> (old_timestep_F); 01743 01744 01745 for (unsigned int i=0; i<dofs_per_cell; ++i) 01746 { 01747 const Tensor<2,dim> pI_LinP = ALETransformations 01748 ::get_pI_LinP<dim> (phi_i_p[i]); 01749 01750 const double J_LinU = ALETransformations 01751 ::get_J_LinU<dim> (q, old_solution_grads, phi_i_grads_u[i]); 01752 01753 const Tensor<2,dim> F_LinU = ALETransformations 01754 ::get_F_LinU<dim> (phi_i_grads_u[i]); 01755 01756 const Tensor<2,dim> F_Inverse_LinU = ALETransformations 01757 ::get_F_Inverse_LinU<dim> 01758 (phi_i_grads_u[i], J, J_LinU, q, old_solution_grads); 01759 01760 const Tensor<2,dim> F_Inverse_T_LinU = transpose(F_Inverse_LinU); 01761 01762 const Tensor<2,dim> J_F_Inverse_T_LinU = ALETransformations 01763 ::get_J_F_Inverse_T_LinU<dim> (phi_i_grads_u[i]); 01764 01765 const Tensor<1,dim> accelaration_term_LinAll = NSEALE 01766 ::get_accelaration_term_LinAll 01767 (phi_i_v[i], v, old_timestep_v, J_LinU, J, old_timestep_J, density_structure); 01768 01769 // STVK: Green-Lagrange strain tensor derivatives 01770 const Tensor<2,dim> E_LinU = 0.5 * (transpose(F_LinU) * F + transpose(F) * F_LinU); 01771 01772 const double tr_E_LinU = StructureTermsALE 01773 ::get_tr_E_LinU<dim> (q,old_solution_grads, phi_i_grads_u[i]); 01774 01775 01776 // STVK 01777 // Piola-kirchhoff stress structure STVK linearized in all directions 01778 Tensor<2,dim> piola_kirchhoff_stress_structure_STVK_LinALL; 01779 piola_kirchhoff_stress_structure_STVK_LinALL = lame_coefficient_lambda * 01780 (F_LinU * tr_E * Identity + F * tr_E_LinU * Identity) 01781 + 2 * lame_coefficient_mu * (F_LinU * E + F * E_LinU); 01782 01783 01784 for (unsigned int j=0; j<dofs_per_cell; ++j) 01785 { 01786 // STVK 01787 const unsigned int comp_j = fe.system_to_component_index(j).first; 01788 if (comp_j == 0 || comp_j == 1) 01789 { 01790 local_matrix(j,i) += (density_structure * phi_i_v[i] * phi_i_v[j] + 01791 timestep * theta * scalar_product(piola_kirchhoff_stress_structure_STVK_LinALL, 01792 phi_i_grads_v[j]) 01793 ) * fe_values.JxW(q); 01794 } 01795 else if (comp_j == 2 || comp_j == 3) 01796 { 01797 local_matrix(j,i) += (density_structure * 1.0/(cell_diameter*cell_diameter) * 01798 (phi_i_u[i] * phi_i_u[j] - timestep * theta * phi_i_v[i] * phi_i_u[j]) 01799 ) * fe_values.JxW(q); 01800 } 01801 else if (comp_j == 4) 01802 { 01803 local_matrix(j,i) += (phi_i_p[i] * phi_i_p[j]) * fe_values.JxW(q); 01804 } 01805 else if (comp_j == 5 || comp_j == 6) 01806 { 01807 local_matrix(j,i) += (alpha_w * (phi_i_w[i] * phi_i_w[j] - scalar_product(phi_i_grads_u[i],phi_i_grads_w[j]) ) 01808 ) * fe_values.JxW(q); 01809 } 01810 // end j dofs 01811 } 01812 // end i dofs 01813 } 01814 // end n_q_points 01815 } 01816 01817 01818 cell->get_dof_indices (local_dof_indices); 01819 constraints.distribute_local_to_global (local_matrix, local_dof_indices, 01820 system_matrix); 01821 // end if (second PDE: STVK material) 01822 } 01823 // end cell 01824 } 01825 01826 timer.exit_section(); 01827 } 01828 01829 01830 01836 template <int dim> 01837 void 01838 FSIALEProblem<dim>::assemble_system_rhs () 01839 { 01840 timer.enter_section("Assemble Rhs."); 01841 system_rhs=0; 01842 01843 QGauss<dim> quadrature_formula(degree+2); 01844 QGauss<dim-1> face_quadrature_formula(degree+2); 01845 01846 FEValues<dim> fe_values (fe, quadrature_formula, 01847 update_values | 01848 update_quadrature_points | 01849 update_JxW_values | 01850 update_gradients); 01851 01852 FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula, 01853 update_values | update_quadrature_points | 01854 update_normal_vectors | update_gradients | 01855 update_JxW_values); 01856 01857 const unsigned int dofs_per_cell = fe.dofs_per_cell; 01858 01859 const unsigned int n_q_points = quadrature_formula.size(); 01860 const unsigned int n_face_q_points = face_quadrature_formula.size(); 01861 01862 Vector<double> local_rhs (dofs_per_cell); 01863 01864 std::vector<unsigned int> local_dof_indices (dofs_per_cell); 01865 01866 const FEValuesExtractors::Vector velocities (0); 01867 const FEValuesExtractors::Vector displacements (dim); 01868 const FEValuesExtractors::Scalar pressure (dim+dim); 01869 const FEValuesExtractors::Vector displacements_w (dim+dim+1); 01870 01871 std::vector<Vector<double> > 01872 old_solution_values (n_q_points, Vector<double>(dim+dim+dim+1)); 01873 01874 std::vector<std::vector<Tensor<1,dim> > > 01875 old_solution_grads (n_q_points, std::vector<Tensor<1,dim> > (dim+dim+dim+1)); 01876 01877 01878 std::vector<Vector<double> > 01879 old_solution_face_values (n_face_q_points, Vector<double>(dim+dim+dim+1)); 01880 01881 std::vector<std::vector<Tensor<1,dim> > > 01882 old_solution_face_grads (n_face_q_points, std::vector<Tensor<1,dim> > (dim+dim+dim+1)); 01883 01884 std::vector<Vector<double> > 01885 old_timestep_solution_values (n_q_points, Vector<double>(dim+dim+dim+1)); 01886 01887 std::vector<std::vector<Tensor<1,dim> > > 01888 old_timestep_solution_grads (n_q_points, std::vector<Tensor<1,dim> > (dim+dim+dim+1)); 01889 01890 std::vector<Vector<double> > 01891 old_timestep_solution_face_values (n_face_q_points, Vector<double>(dim+dim+dim+1)); 01892 01893 std::vector<std::vector<Tensor<1,dim> > > 01894 old_timestep_solution_face_grads (n_face_q_points, std::vector<Tensor<1,dim> > (dim+dim+dim+1)); 01895 01896 01897 typename DoFHandler<dim>::active_cell_iterator 01898 cell = dof_handler.begin_active(), 01899 endc = dof_handler.end(); 01900 01901 for (; cell!=endc; ++cell) 01902 { 01903 fe_values.reinit (cell); 01904 local_rhs = 0; 01905 01906 cell_diameter = cell->diameter(); 01907 01908 // old Newton iteration 01909 fe_values.get_function_values (solution, old_solution_values); 01910 fe_values.get_function_grads (solution, old_solution_grads); 01911 01912 // old timestep iteration 01913 fe_values.get_function_values (old_timestep_solution, old_timestep_solution_values); 01914 fe_values.get_function_grads (old_timestep_solution, old_timestep_solution_grads); 01915 01916 // Again, material_id == 0 corresponds to 01917 // the domain for fluid equations 01918 if (cell->material_id() == 0) 01919 { 01920 for (unsigned int q=0; q<n_q_points; ++q) 01921 { 01922 const Tensor<2,dim> pI = ALETransformations 01923 ::get_pI<dim> (q, old_solution_values); 01924 01925 const Tensor<1,dim> v = ALETransformations 01926 ::get_v<dim> (q, old_solution_values); 01927 01928 const Tensor<2,dim> grad_v = ALETransformations 01929 ::get_grad_v<dim> (q, old_solution_grads); 01930 01931 const Tensor<2,dim> grad_u = ALETransformations 01932 ::get_grad_u<dim> (q, old_solution_grads); 01933 01934 const Tensor<2,dim> grad_v_T = ALETransformations 01935 ::get_grad_v_T<dim> (grad_v); 01936 01937 const Tensor<1,dim> u = ALETransformations 01938 ::get_u<dim> (q, old_solution_values); 01939 01940 const Tensor<1,dim> w = ALETransformations 01941 ::get_w<dim> (q, old_solution_values); 01942 01943 const Tensor<2,dim> grad_w = ALETransformations 01944 ::get_grad_w<dim> (q, old_solution_grads); 01945 01946 const Tensor<2,dim> F = ALETransformations 01947 ::get_F<dim> (q, old_solution_grads); 01948 01949 const Tensor<2,dim> F_Inverse = ALETransformations 01950 ::get_F_Inverse<dim> (F); 01951 01952 const Tensor<2,dim> F_Inverse_T = ALETransformations 01953 ::get_F_Inverse_T<dim> (F_Inverse); 01954 01955 const double J = ALETransformations 01956 ::get_J<dim> (F); 01957 01958 // This is the fluid stress tensor in ALE formulation 01959 const Tensor<2,dim> sigma_ALE = NSEALE 01960 ::get_stress_fluid_except_pressure_ALE<dim> 01961 (density_fluid, viscosity, grad_v, grad_v_T, F_Inverse, F_Inverse_T ); 01962 01963 // We proceed by catching the previous time step values 01964 const Tensor<2,dim> old_timestep_pI = ALETransformations 01965 ::get_pI<dim> (q, old_timestep_solution_values); 01966 01967 const Tensor<1,dim> old_timestep_v = ALETransformations 01968 ::get_v<dim> (q, old_timestep_solution_values); 01969 01970 const Tensor<2,dim> old_timestep_grad_v = ALETransformations 01971 ::get_grad_v<dim> (q, old_timestep_solution_grads); 01972 01973 const Tensor<2,dim> old_timestep_grad_v_T = ALETransformations 01974 ::get_grad_v_T<dim> (old_timestep_grad_v); 01975 01976 const Tensor<1,dim> old_timestep_u = ALETransformations 01977 ::get_u<dim> (q, old_timestep_solution_values); 01978 01979 const Tensor<2,dim> old_timestep_grad_u = ALETransformations 01980 ::get_grad_u<dim> (q, old_timestep_solution_grads); 01981 01982 const Tensor<2,dim> old_timestep_F = ALETransformations 01983 ::get_F<dim> (q, old_timestep_solution_grads); 01984 01985 const Tensor<2,dim> old_timestep_F_Inverse = ALETransformations 01986 ::get_F_Inverse<dim> (old_timestep_F); 01987 01988 const Tensor<2,dim> old_timestep_F_Inverse_T = ALETransformations 01989 ::get_F_Inverse_T<dim> (old_timestep_F_Inverse); 01990 01991 const double old_timestep_J = ALETransformations 01992 ::get_J<dim> (old_timestep_F); 01993 01994 // This is the fluid stress tensor in the ALE formulation 01995 // at the previous time step 01996 const Tensor<2,dim> old_timestep_sigma_ALE = NSEALE 01997 ::get_stress_fluid_except_pressure_ALE<dim> 01998 (density_fluid, viscosity, old_timestep_grad_v, old_timestep_grad_v_T, 01999 old_timestep_F_Inverse, old_timestep_F_Inverse_T ); 02000 02001 Tensor<2,dim> stress_fluid; 02002 stress_fluid.clear(); 02003 stress_fluid = (J * sigma_ALE * F_Inverse_T); 02004 02005 Tensor<2,dim> fluid_pressure; 02006 fluid_pressure.clear(); 02007 fluid_pressure = (-pI * J * F_Inverse_T); 02008 02009 Tensor<2,dim> old_timestep_stress_fluid; 02010 old_timestep_stress_fluid.clear(); 02011 old_timestep_stress_fluid = 02012 (old_timestep_J * old_timestep_sigma_ALE * old_timestep_F_Inverse_T); 02013 02014 // Divergence of the fluid in the ALE formulation 02015 const double incompressiblity_fluid = NSEALE 02016 ::get_Incompressibility_ALE<dim> (q, old_solution_grads); 02017 02018 // Convection term of the fluid in the ALE formulation. 02019 // We emphasize that the fluid convection term for 02020 // non-stationary flow problems in ALE 02021 // representation is difficult to derive. 02022 // For adequate discretization, the convection term will be 02023 // split into three smaller terms: 02024 Tensor<1,dim> convection_fluid; 02025 convection_fluid.clear(); 02026 convection_fluid = density_fluid * J * (grad_v * F_Inverse * v); 02027 02028 // The second convection term for the fluid in the ALE formulation 02029 Tensor<1,dim> convection_fluid_with_u; 02030 convection_fluid_with_u.clear(); 02031 convection_fluid_with_u = 02032 density_fluid * J * (grad_v * F_Inverse * u); 02033 02034 // The third convection term for the fluid in the ALE formulation 02035 Tensor<1,dim> convection_fluid_with_old_timestep_u; 02036 convection_fluid_with_old_timestep_u.clear(); 02037 convection_fluid_with_old_timestep_u = 02038 density_fluid * J * (grad_v * F_Inverse * old_timestep_u); 02039 02040 // The convection term of the previous time step 02041 Tensor<1,dim> old_timestep_convection_fluid; 02042 old_timestep_convection_fluid.clear(); 02043 old_timestep_convection_fluid = 02044 (density_fluid * old_timestep_J * 02045 (old_timestep_grad_v * old_timestep_F_Inverse * old_timestep_v)); 02046 02047 for (unsigned int i=0; i<dofs_per_cell; ++i) 02048 { 02049 // Fluid, NSE in ALE 02050 const unsigned int comp_i = fe.system_to_component_index(i).first; 02051 if (comp_i == 0 || comp_i == 1) 02052 { 02053 const Tensor<1,dim> phi_i_v = fe_values[velocities].value (i, q); 02054 const Tensor<2,dim> phi_i_grads_v = fe_values[velocities].gradient (i, q); 02055 02056 local_rhs(i) -= (density_fluid * (J + old_timestep_J)/2.0 * 02057 (v - old_timestep_v) * phi_i_v + 02058 timestep * theta * convection_fluid * phi_i_v + 02059 timestep * (1.0-theta) * 02060 old_timestep_convection_fluid * phi_i_v - 02061 (convection_fluid_with_u - 02062 convection_fluid_with_old_timestep_u) * phi_i_v + 02063 timestep * scalar_product(fluid_pressure, phi_i_grads_v) + 02064 timestep * theta * scalar_product(stress_fluid, phi_i_grads_v) + 02065 timestep * (1.0-theta) * 02066 scalar_product(old_timestep_stress_fluid, phi_i_grads_v) 02067 ) * fe_values.JxW(q); 02068 02069 } 02070 else if (comp_i == 2 || comp_i == 3) 02071 { 02072 const Tensor<1,dim> phi_i_u = fe_values[displacements].value (i, q); 02073 const Tensor<2,dim> phi_i_grads_u = fe_values[displacements].gradient (i, q); 02074 02075 local_rhs(i) -= (alpha_u * scalar_product(grad_w, phi_i_grads_u) 02076 ) * fe_values.JxW(q); 02077 } 02078 else if (comp_i == 4) 02079 { 02080 const double phi_i_p = fe_values[pressure].value (i, q); 02081 local_rhs(i) -= (incompressiblity_fluid * phi_i_p) * fe_values.JxW(q); 02082 } 02083 else if (comp_i == 5 || comp_i == 6) 02084 { 02085 const Tensor<1,dim> phi_i_w = fe_values[displacements_w].value (i, q); 02086 const Tensor<2,dim> phi_i_grads_w = fe_values[displacements_w].gradient (i, q); 02087 02088 local_rhs(i) -= alpha_w * (w * phi_i_w - scalar_product(grad_u,phi_i_grads_w)) * 02089 fe_values.JxW(q); 02090 } 02091 // end i dofs 02092 } 02093 // close n_q_points 02094 } 02095 02096 // As already discussed in the assembling method for the matrix, 02097 // we have to integrate some terms on the outflow boundary: 02098 for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face) 02099 { 02100 if (cell->face(face)->at_boundary() && 02101 (cell->face(face)->boundary_indicator() == 1) 02102 ) 02103 { 02104 02105 fe_face_values.reinit (cell, face); 02106 02107 fe_face_values.get_function_values (solution, old_solution_face_values); 02108 fe_face_values.get_function_grads (solution, old_solution_face_grads); 02109 02110 fe_face_values.get_function_values (old_timestep_solution, old_timestep_solution_face_values); 02111 fe_face_values.get_function_grads (old_timestep_solution, old_timestep_solution_face_grads); 02112 02113 for (unsigned int q=0; q<n_face_q_points; ++q) 02114 { 02115 // These are terms coming from the 02116 // previous Newton iterations ... 02117 const Tensor<1,dim> v = ALETransformations 02118 ::get_v<dim> (q, old_solution_face_values); 02119 02120 const Tensor<2,dim> grad_v = ALETransformations 02121 ::get_grad_v<dim> (q, old_solution_face_grads); 02122 02123 const Tensor<2,dim> grad_v_T = ALETransformations 02124 ::get_grad_v_T<dim> (grad_v); 02125 02126 const Tensor<2,dim> F = ALETransformations 02127 ::get_F<dim> (q, old_solution_face_grads); 02128 02129 const Tensor<2,dim> F_Inverse = ALETransformations 02130 ::get_F_Inverse<dim> (F); 02131 02132 const Tensor<2,dim> F_Inverse_T = ALETransformations 02133 ::get_F_Inverse_T<dim> (F_Inverse); 02134 02135 const double J = ALETransformations 02136 ::get_J<dim> (F); 02137 02138 // ... and here from the previous time step iteration 02139 const Tensor<1,dim> old_timestep_v = ALETransformations 02140 ::get_v<dim> (q, old_timestep_solution_face_values); 02141 02142 const Tensor<2,dim> old_timestep_grad_v = ALETransformations 02143 ::get_grad_v<dim> (q, old_timestep_solution_face_grads); 02144 02145 const Tensor<2,dim> old_timestep_grad_v_T = ALETransformations 02146 ::get_grad_v_T<dim> (old_timestep_grad_v); 02147 02148 const Tensor<2,dim> old_timestep_F = ALETransformations 02149 ::get_F<dim> (q, old_timestep_solution_face_grads); 02150 02151 const Tensor<2,dim> old_timestep_F_Inverse = ALETransformations 02152 ::get_F_Inverse<dim> (old_timestep_F); 02153 02154 const Tensor<2,dim> old_timestep_F_Inverse_T = ALETransformations 02155 ::get_F_Inverse_T<dim> (old_timestep_F_Inverse); 02156 02157 const double old_timestep_J = ALETransformations 02158 ::get_J<dim> (old_timestep_F); 02159 02160 Tensor<2,dim> sigma_ALE_tilde; 02161 sigma_ALE_tilde.clear(); 02162 sigma_ALE_tilde = 02163 (density_fluid * viscosity * F_Inverse_T * grad_v_T); 02164 02165 Tensor<2,dim> old_timestep_sigma_ALE_tilde; 02166 old_timestep_sigma_ALE_tilde.clear(); 02167 old_timestep_sigma_ALE_tilde = 02168 (density_fluid * viscosity * old_timestep_F_Inverse_T * old_timestep_grad_v_T); 02169 02170 // Neumann boundary integral 02171 Tensor<2,dim> stress_fluid_transposed_part; 02172 stress_fluid_transposed_part.clear(); 02173 stress_fluid_transposed_part = (J * sigma_ALE_tilde * F_Inverse_T); 02174 02175 Tensor<2,dim> old_timestep_stress_fluid_transposed_part; 02176 old_timestep_stress_fluid_transposed_part.clear(); 02177 old_timestep_stress_fluid_transposed_part = 02178 (old_timestep_J * old_timestep_sigma_ALE_tilde * old_timestep_F_Inverse_T); 02179 02180 const Tensor<1,dim> neumann_value 02181 = (stress_fluid_transposed_part * fe_face_values.normal_vector(q)); 02182 02183 const Tensor<1,dim> old_timestep_neumann_value 02184 = (old_timestep_stress_fluid_transposed_part * fe_face_values.normal_vector(q)); 02185 02186 for (unsigned int i=0; i<dofs_per_cell; ++i) 02187 { 02188 const unsigned int comp_i = fe.system_to_component_index(i).first; 02189 if (comp_i == 0 || comp_i == 1) 02190 { 02191 local_rhs(i) += (timestep * theta * 02192 neumann_value * fe_face_values[velocities].value (i, q) + 02193 timestep * (1.0-theta) * 02194 old_timestep_neumann_value * 02195 fe_face_values[velocities].value (i, q) 02196 ) * fe_face_values.JxW(q); 02197 } 02198 // end i 02199 } 02200 // end face_n_q_points 02201 } 02202 } 02203 } // end face integrals do-nothing condition 02204 02205 02206 // The computation of these face integrals on the interface has 02207 // already been discussed in the matrix section. 02208 for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face) 02209 { 02210 if (cell->neighbor_index(face) != -1) 02211 if (cell->material_id() != cell->neighbor(face)->material_id() && 02212 cell->face(face)->boundary_indicator()!=81) 02213 { 02214 fe_face_values.reinit (cell, face); 02215 fe_face_values.get_function_grads (solution, old_solution_face_grads); 02216 02217 for (unsigned int q=0; q<n_face_q_points; ++q) 02218 { 02219 const Tensor<2,dim> grad_u = ALETransformations 02220 ::get_grad_u<dim> (q, old_solution_face_grads); 02221 02222 const Tensor<2,dim> grad_w = ALETransformations 02223 ::get_grad_w<dim> (q, old_solution_face_grads); 02224 02225 const Tensor<1,dim> neumann_value_u 02226 = (grad_u * fe_face_values.normal_vector(q)); 02227 02228 const Tensor<1,dim> neumann_value_w 02229 = (grad_w * fe_face_values.normal_vector(q)); 02230 02231 for (unsigned int i=0; i<dofs_per_cell; ++i) 02232 { 02233 const unsigned int comp_i = fe.system_to_component_index(i).first; 02234 if (comp_i == 2 || comp_i == 3) 02235 { 02236 local_rhs(i) += (alpha_w * neumann_value_w * 02237 fe_face_values[displacements].value (i, q)) * 02238 fe_face_values.JxW(q); 02239 } 02240 else if (comp_i == 5 || comp_i == 6) 02241 { 02242 local_rhs(i) += (alpha_u * neumann_value_u * 02243 fe_face_values[displacements_w].value (i, q)) * 02244 fe_face_values.JxW(q); 02245 } 02246 02247 } // end i 02248 02249 } // end face_n_q_points 02250 02251 } 02252 } // end face for interface conditions 02253 02254 02255 cell->get_dof_indices (local_dof_indices); 02256 constraints.distribute_local_to_global (local_rhs, local_dof_indices, 02257 system_rhs); 02258 02259 // Finally, we arrive at the end for assembling 02260 // the variational formulation for the fluid part and step to 02261 // the assembling process of the structure terms: 02262 } 02263 else if (cell->material_id() == 1) 02264 { 02265 for (unsigned int q=0; q<n_q_points; ++q) 02266 { 02267 const Tensor<2,dim> pI = ALETransformations 02268 ::get_pI<dim> (q, old_solution_values); 02269 02270 const Tensor<1,dim> v = ALETransformations 02271 ::get_v<dim> (q, old_solution_values); 02272 02273 const Tensor<2,dim> grad_v = ALETransformations 02274 ::get_grad_v<dim> (q, old_solution_grads); 02275 02276 const Tensor<2,dim> grad_v_T = ALETransformations 02277 ::get_grad_v_T<dim> (grad_v); 02278 02279 const Tensor<1,dim> u = ALETransformations 02280 ::get_u<dim> (q, old_solution_values); 02281 02282 const Tensor<2,dim> grad_u = ALETransformations 02283 ::get_grad_u<dim> (q, old_solution_grads); 02284 02285 const Tensor<1,dim> w = ALETransformations 02286 ::get_w<dim> (q, old_solution_values); 02287 02288 const Tensor<2,dim> grad_w = ALETransformations 02289 ::get_grad_w<dim> (q, old_solution_grads); 02290 02291 const Tensor<2,dim> F = ALETransformations 02292 ::get_F<dim> (q, old_solution_grads); 02293 02294 const Tensor<2,dim> F_T = ALETransformations 02295 ::get_F_T<dim> (F); 02296 02297 const Tensor<2,dim> Identity = ALETransformations 02298 ::get_Identity<dim> (); 02299 02300 const Tensor<2,dim> F_Inverse = ALETransformations 02301 ::get_F_Inverse<dim> (F); 02302 02303 const Tensor<2,dim> F_Inverse_T = ALETransformations 02304 ::get_F_Inverse_T<dim> (F_Inverse); 02305 02306 const double J = ALETransformations 02307 ::get_J<dim> (F); 02308 02309 const Tensor<2,dim> E = StructureTermsALE 02310 ::get_E<dim> (F_T, F, Identity); 02311 02312 const double tr_E = StructureTermsALE 02313 ::get_tr_E<dim> (E); 02314 02315 // Previous time step values 02316 const Tensor<2,dim> old_timestep_pI = ALETransformations 02317 ::get_pI<dim> (q, old_timestep_solution_values); 02318 02319 const Tensor<1,dim> old_timestep_v = ALETransformations 02320 ::get_v<dim> (q, old_timestep_solution_values); 02321 02322 const Tensor<2,dim> old_timestep_grad_v = ALETransformations 02323 ::get_grad_v<dim> (q, old_timestep_solution_grads); 02324 02325 const Tensor<2,dim> old_timestep_grad_v_T = ALETransformations 02326 ::get_grad_v_T<dim> (old_timestep_grad_v); 02327 02328 const Tensor<1,dim> old_timestep_u = ALETransformations 02329 ::get_u<dim> (q, old_timestep_solution_values); 02330 02331 const Tensor<2,dim> old_timestep_F = ALETransformations 02332 ::get_F<dim> (q, old_timestep_solution_grads); 02333 02334 const Tensor<2,dim> old_timestep_F_Inverse = ALETransformations 02335 ::get_F_Inverse<dim> (old_timestep_F); 02336 02337 const Tensor<2,dim> old_timestep_F_T = ALETransformations 02338 ::get_F_T<dim> (old_timestep_F); 02339 02340 const Tensor<2,dim> old_timestep_F_Inverse_T = ALETransformations 02341 ::get_F_Inverse_T<dim> (old_timestep_F_Inverse); 02342 02343 const double old_timestep_J = ALETransformations 02344 ::get_J<dim> (old_timestep_F); 02345 02346 const Tensor<2,dim> old_timestep_E = StructureTermsALE 02347 ::get_E<dim> (old_timestep_F_T, old_timestep_F, Identity); 02348 02349 const double old_timestep_tr_E = StructureTermsALE 02350 ::get_tr_E<dim> (old_timestep_E); 02351 02352 02353 // STVK structure model 02354 Tensor<2,dim> sigma_structure_ALE; 02355 sigma_structure_ALE.clear(); 02356 sigma_structure_ALE = (1.0/J * 02357 F * (lame_coefficient_lambda * 02358 tr_E * Identity + 02359 2 * lame_coefficient_mu * 02360 E) * 02361 F_T); 02362 02363 02364 Tensor<2,dim> stress_term; 02365 stress_term.clear(); 02366 stress_term = (J * sigma_structure_ALE * F_Inverse_T); 02367 02368 Tensor<2,dim> old_timestep_sigma_structure_ALE; 02369 old_timestep_sigma_structure_ALE.clear(); 02370 old_timestep_sigma_structure_ALE = (1.0/old_timestep_J * 02371 old_timestep_F * (lame_coefficient_lambda * 02372 old_timestep_tr_E * Identity + 02373 2 * lame_coefficient_mu * 02374 old_timestep_E) * 02375 old_timestep_F_T); 02376 02377 Tensor<2,dim> old_timestep_stress_term; 02378 old_timestep_stress_term.clear(); 02379 old_timestep_stress_term = (old_timestep_J * old_timestep_sigma_structure_ALE * old_timestep_F_Inverse_T); 02380 02381 for (unsigned int i=0; i<dofs_per_cell; ++i) 02382 { 02383 // STVK structure model 02384 const unsigned int comp_i = fe.system_to_component_index(i).first; 02385 if (comp_i == 0 || comp_i == 1) 02386 { 02387 const Tensor<1,dim> phi_i_v = fe_values[velocities].value (i, q); 02388 const Tensor<2,dim> phi_i_grads_v = fe_values[velocities].gradient (i, q); 02389 02390 local_rhs(i) -= (density_structure * (v - old_timestep_v) * phi_i_v + 02391 timestep * theta * scalar_product(stress_term,phi_i_grads_v) + 02392 timestep * (1.0-theta) * scalar_product(old_timestep_stress_term, phi_i_grads_v) 02393 ) * fe_values.JxW(q); 02394 02395 } 02396 else if (comp_i == 2 || comp_i == 3) 02397 { 02398 const Tensor<1,dim> phi_i_u = fe_values[displacements].value (i, q); 02399 local_rhs(i) -= (density_structure * 1.0/(cell_diameter*cell_diameter) * 02400 ((u - old_timestep_u) * phi_i_u - 02401 timestep * (theta * v + (1.0-theta) * 02402 old_timestep_v) * phi_i_u) 02403 ) * fe_values.JxW(q); 02404 02405 } 02406 else if (comp_i == 4) 02407 { 02408 const double phi_i_p = fe_values[pressure].value (i, q); 02409 local_rhs(i) -= (old_solution_values[q](dim+dim) * phi_i_p) * fe_values.JxW(q); 02410 02411 } 02412 else if (comp_i == 5 || comp_i == 6) 02413 { 02414 const Tensor<1,dim> phi_i_w = fe_values[displacements_w].value (i, q); 02415 const Tensor<2,dim> phi_i_grads_w = fe_values[displacements_w].gradient (i, q); 02416 02417 local_rhs(i) -= alpha_w * (w * phi_i_w - scalar_product(grad_u,phi_i_grads_w)) * 02418 fe_values.JxW(q); 02419 } 02420 // end i 02421 } 02422 // end n_q_points 02423 } 02424 02425 cell->get_dof_indices (local_dof_indices); 02426 constraints.distribute_local_to_global (local_rhs, local_dof_indices, 02427 system_rhs); 02428 02429 // end if (for STVK material) 02430 } 02431 02432 } // end cell 02433 02434 timer.exit_section(); 02435 } 02436 02437 02455 template <int dim> 02456 void 02457 FSIALEProblem<dim>::set_initial_bc (const double time) 02458 { 02459 std::map<unsigned int,double> boundary_values; 02460 std::vector<bool> component_mask (dim+dim+dim+1, true); 02461 // (Scalar) pressure 02462 component_mask[dim+dim] = false; 02463 02464 // Additional displacement w: 02465 component_mask[dim+dim+1] = false; 02466 component_mask[dim+dim+dim] = false; 02467 02468 VectorTools::interpolate_boundary_values (dof_handler, 02469 0, 02470 BoundaryParabolic<dim>(time), 02471 boundary_values, 02472 component_mask); 02473 02474 VectorTools::interpolate_boundary_values (dof_handler, 02475 2, 02476 ZeroFunction<dim>(dim+dim+dim+1), 02477 boundary_values, 02478 component_mask); 02479 02480 VectorTools::interpolate_boundary_values (dof_handler, 02481 80, 02482 ZeroFunction<dim>(dim+dim+dim+1), 02483 boundary_values, 02484 component_mask); 02485 02486 VectorTools::interpolate_boundary_values (dof_handler, 02487 81, 02488 ZeroFunction<dim>(dim+dim+dim+1), 02489 boundary_values, 02490 component_mask); 02491 02492 component_mask[0] = false; 02493 component_mask[1] = false; 02494 02495 VectorTools::interpolate_boundary_values (dof_handler, 02496 1, 02497 ZeroFunction<dim>(dim+dim+dim+1), 02498 boundary_values, 02499 component_mask); 02500 02501 for (typename std::map<unsigned int, double>::const_iterator 02502 i = boundary_values.begin(); 02503 i != boundary_values.end(); 02504 ++i) 02505 solution(i->first) = i->second; 02506 02507 } 02508 02515 template <int dim> 02516 void 02517 FSIALEProblem<dim>::set_newton_bc () 02518 { 02519 std::vector<bool> component_mask (dim+dim+dim+1, true); 02520 component_mask[dim+dim] = false; 02521 component_mask[dim+dim+1] = false; 02522 component_mask[dim+dim+dim] = false; 02523 02524 VectorTools::interpolate_boundary_values (dof_handler, 02525 0, 02526 ZeroFunction<dim>(dim+dim+dim+1), 02527 constraints, 02528 component_mask); 02529 02530 VectorTools::interpolate_boundary_values (dof_handler, 02531 2, 02532 ZeroFunction<dim>(dim+dim+dim+1), 02533 constraints, 02534 component_mask); 02535 02536 VectorTools::interpolate_boundary_values (dof_handler, 02537 80, 02538 ZeroFunction<dim>(dim+dim+dim+1), 02539 constraints, 02540 component_mask); 02541 VectorTools::interpolate_boundary_values (dof_handler, 02542 81, 02543 ZeroFunction<dim>(dim+dim+dim+1), 02544 constraints, 02545 component_mask); 02546 component_mask[0] = false; 02547 component_mask[1] = false; 02548 02549 VectorTools::interpolate_boundary_values (dof_handler, 02550 1, 02551 ZeroFunction<dim>(dim+dim+dim+1), 02552 constraints, 02553 component_mask); 02554 } 02555 02556 02562 template <int dim> 02563 void FSIALEProblem<dim>::newton_iteration (const double time) 02564 02565 { 02566 Timer timer_newton; 02567 const double lower_bound_newton_residuum = 1.0e-8; 02568 const unsigned int max_no_newton_steps = 40; 02569 02570 // Decision whether the system matrix should be build 02571 // at each Newton step 02572 const double nonlinear_rho = 0.1; 02573 02574 // Line search parameters 02575 unsigned int line_search_step; 02576 const unsigned int max_no_line_search_steps = 10; 02577 const double line_search_damping = 0.6; 02578 double new_newton_residuum; 02579 02580 // Application of the initial boundary conditions to the 02581 // variational equations: 02582 set_initial_bc (time); 02583 assemble_system_rhs(); 02584 02585 double newton_residuum = system_rhs.linfty_norm(); 02586 double old_newton_residuum= newton_residuum; 02587 unsigned int newton_step = 1; 02588 02589 if (newton_residuum < lower_bound_newton_residuum) 02590 { 02591 std::cout << '\t' 02592 << std::scientific 02593 << newton_residuum 02594 << std::endl; 02595 } 02596 02597 while (newton_residuum > lower_bound_newton_residuum && 02598 newton_step < max_no_newton_steps) 02599 { 02600 timer_newton.start(); 02601 old_newton_residuum = newton_residuum; 02602 02603 assemble_system_rhs(); 02604 newton_residuum = system_rhs.linfty_norm(); 02605 02606 if (newton_residuum < lower_bound_newton_residuum) 02607 { 02608 std::cout << '\t' 02609 << std::scientific 02610 << newton_residuum << std::endl; 02611 break; 02612 } 02613 02614 if (newton_residuum/old_newton_residuum > nonlinear_rho) 02615 assemble_system_matrix (); 02616 02617 // Solve Ax = b 02618 solve (); 02619 02620 line_search_step = 0; 02621 for ( ; 02622 line_search_step < max_no_line_search_steps; 02623 ++line_search_step) 02624 { 02625 solution += newton_update; 02626 02627 assemble_system_rhs (); 02628 new_newton_residuum = system_rhs.linfty_norm(); 02629 02630 if (new_newton_residuum < newton_residuum) 02631 break; 02632 else 02633 solution -= newton_update; 02634 02635 newton_update *= line_search_damping; 02636 } 02637 02638 timer_newton.stop(); 02639 02640 std::cout << std::setprecision(5) <<newton_step << '\t' 02641 << std::scientific << newton_residuum << '\t' 02642 << std::scientific << newton_residuum/old_newton_residuum <<'\t' ; 02643 if (newton_residuum/old_newton_residuum > nonlinear_rho) 02644 std::cout << "r" << '\t' ; 02645 else 02646 std::cout << " " << '\t' ; 02647 std::cout << line_search_step << '\t' 02648 << std::scientific << timer_newton () 02649 << std::endl; 02650 02651 02652 // Updates 02653 timer_newton.reset(); 02654 newton_step++; 02655 } 02656 } 02657 02676 template <int dim> 02677 void 02678 FSIALEProblem<dim>::solve () 02679 { 02680 timer.enter_section("Solve linear system."); 02681 Vector<double> sol, rhs; 02682 sol = newton_update; 02683 rhs = system_rhs; 02684 02685 SparseDirectUMFPACK A_direct; 02686 A_direct.factorize(system_matrix); 02687 A_direct.vmult(sol,rhs); 02688 newton_update = sol; 02689 02690 constraints.distribute (newton_update); 02691 timer.exit_section(); 02692 } 02693 02694 02695 /* This function is known from almost all other 02696 * tutorial steps in deal.II. 02697 * However, we emphasize that the FSI problem 02698 * is computed on a fixed mesh (instead of 02699 * moving the mesh as done in other references; 02700 * we refer the reader to the accompanying article 02701 * and the comments made therein). 02702 * For this reason, the output of the 02703 * solution in *.vtk format corresponds to the 02704 * solution on the fixed mesh and, therefore, 02705 * in a visualization program, the reader 02706 * has to postprocess the solution. 02707 * We also refer to the MappingQEulerian in 02708 * the deal.II librar, which is able 02709 * to tansform the solution to the current (i.e., 02710 * the physical mesh) 02711 */ 02712 template <int dim> 02713 void 02714 FSIALEProblem<dim>::output_results (const unsigned int refinement_cycle, 02715 const BlockVector<double> output_vector) const 02716 { 02717 02718 std::vector<std::string> solution_names (dim, "velocity"); 02719 solution_names.push_back ("displacement"); 02720 solution_names.push_back ("displacement"); 02721 solution_names.push_back ("p_fluid"); 02722 solution_names.push_back ("displace_w"); 02723 solution_names.push_back ("displace_w"); 02724 02725 std::vector<DataComponentInterpretation::DataComponentInterpretation> 02726 data_component_interpretation 02727 (dim+dim, DataComponentInterpretation::component_is_part_of_vector); 02728 02729 data_component_interpretation 02730 .push_back (DataComponentInterpretation::component_is_scalar); 02731 02732 data_component_interpretation 02733 .push_back (DataComponentInterpretation::component_is_part_of_vector); 02734 data_component_interpretation 02735 .push_back (DataComponentInterpretation::component_is_part_of_vector); 02736 02737 DataOut<dim> data_out; 02738 data_out.attach_dof_handler (dof_handler); 02739 02740 data_out.add_data_vector (output_vector, solution_names, 02741 DataOut<dim>::type_dof_data, 02742 data_component_interpretation); 02743 02744 data_out.build_patches (); 02745 02746 std::string filename_basis; 02747 filename_basis = "solution_fsi_1_"; 02748 02749 std::ostringstream filename; 02750 02751 std::cout << "------------------" << std::endl; 02752 std::cout << "Write solution" << std::endl; 02753 std::cout << "------------------" << std::endl; 02754 std::cout << std::endl; 02755 filename << filename_basis 02756 << Utilities::int_to_string (refinement_cycle, 5) 02757 << ".vtk"; 02758 02759 std::ofstream output (filename.str().c_str()); 02760 data_out.write_vtk (output); 02761 02762 } 02763 02769 template <int dim> 02770 double FSIALEProblem<dim>::compute_point_value (Point<dim> p, 02771 const unsigned int component) const 02772 { 02773 02774 Vector<double> tmp_vector(dim+dim+dim+1); 02775 VectorTools::point_value (dof_handler, 02776 solution, 02777 p, 02778 tmp_vector); 02779 02780 return tmp_vector(component); 02781 } 02782 02790 template <int dim> 02791 void FSIALEProblem<dim>::compute_drag_lift_fsi_fluid_tensor() 02792 { 02793 02794 const QGauss<dim-1> face_quadrature_formula (3); 02795 FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula, 02796 update_values | update_gradients | update_normal_vectors | 02797 update_JxW_values); 02798 02799 const unsigned int dofs_per_cell = fe.dofs_per_cell; 02800 const unsigned int n_face_q_points = face_quadrature_formula.size(); 02801 02802 std::vector<unsigned int> local_dof_indices (dofs_per_cell); 02803 std::vector<Vector<double> > face_solution_values (n_face_q_points, 02804 Vector<double> (dim+dim+dim+1)); 02805 02806 std::vector<std::vector<Tensor<1,dim> > > 02807 face_solution_grads (n_face_q_points, std::vector<Tensor<1,dim> > (dim+dim+dim+1)); 02808 02809 Tensor<1,dim> drag_lift_value; 02810 02811 typename DoFHandler<dim>::active_cell_iterator 02812 cell = dof_handler.begin_active(), 02813 endc = dof_handler.end(); 02814 02815 for (; cell!=endc; ++cell) 02816 { 02817 02818 // First, we are going to compute the forces that 02819 // act on the cylinder. We notice that only the fluid 02820 // equations are defined here. 02821 for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face) 02822 if (cell->face(face)->at_boundary() && 02823 cell->face(face)->boundary_indicator()==80) 02824 { 02825 fe_face_values.reinit (cell, face); 02826 fe_face_values.get_function_values (solution, face_solution_values); 02827 fe_face_values.get_function_grads (solution, face_solution_grads); 02828 02829 for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point) 02830 { 02831 const Tensor<2,dim> pI = ALETransformations 02832 ::get_pI<dim> (q_point, face_solution_values); 02833 02834 const Tensor<1,dim> v = ALETransformations 02835 ::get_v<dim> (q_point, face_solution_values); 02836 02837 const Tensor<2,dim> grad_v = ALETransformations 02838 ::get_grad_v<dim> (q_point, face_solution_grads); 02839 02840 const Tensor<2,dim> grad_v_T = ALETransformations 02841 ::get_grad_v_T<dim> (grad_v); 02842 02843 const Tensor<2,dim> F = ALETransformations 02844 ::get_F<dim> (q_point, face_solution_grads); 02845 02846 const Tensor<2,dim> F_Inverse = ALETransformations 02847 ::get_F_Inverse<dim> (F); 02848 02849 const Tensor<2,dim> F_Inverse_T = ALETransformations 02850 ::get_F_Inverse_T<dim> (F_Inverse); 02851 02852 const double J = ALETransformations 02853 ::get_J<dim> (F); 02854 02855 const Tensor<2,dim> sigma_ALE = NSEALE 02856 ::get_stress_fluid_except_pressure_ALE<dim> 02857 (density_fluid, viscosity, 02858 grad_v, grad_v_T, F_Inverse, F_Inverse_T ); 02859 02860 Tensor<2,dim> stress_fluid; 02861 stress_fluid.clear(); 02862 stress_fluid = (J * sigma_ALE * F_Inverse_T); 02863 02864 Tensor<2,dim> fluid_pressure; 02865 fluid_pressure.clear(); 02866 fluid_pressure = (-pI * J * F_Inverse_T); 02867 02868 drag_lift_value -= (stress_fluid + fluid_pressure) * 02869 fe_face_values.normal_vector(q_point)* fe_face_values.JxW(q_point); 02870 02871 } 02872 } // end boundary 80 for fluid 02873 02874 // Now, we compute the forces that act on the beam. Here, 02875 // we have two possibilities as already discussed in the paper. 02876 // We use again the fluid tensor to compute 02877 // drag and lift: 02878 if (cell->material_id() == 0) 02879 { 02880 for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face) 02881 if (cell->neighbor_index(face) != -1) 02882 if (cell->material_id() != cell->neighbor(face)->material_id() && 02883 cell->face(face)->boundary_indicator()!=80) 02884 { 02885 02886 fe_face_values.reinit (cell, face); 02887 fe_face_values.get_function_values (solution, face_solution_values); 02888 fe_face_values.get_function_grads (solution, face_solution_grads); 02889 02890 for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point) 02891 { 02892 const Tensor<2,dim> pI = ALETransformations 02893 ::get_pI<dim> (q_point, face_solution_values); 02894 02895 const Tensor<1,dim> v = ALETransformations 02896 ::get_v<dim> (q_point, face_solution_values); 02897 02898 const Tensor<2,dim> grad_v = ALETransformations 02899 ::get_grad_v<dim> (q_point, face_solution_grads); 02900 02901 const Tensor<2,dim> grad_v_T = ALETransformations 02902 ::get_grad_v_T<dim> (grad_v); 02903 02904 const Tensor<2,dim> F = ALETransformations 02905 ::get_F<dim> (q_point, face_solution_grads); 02906 02907 const Tensor<2,dim> F_Inverse = ALETransformations 02908 ::get_F_Inverse<dim> (F); 02909 02910 const Tensor<2,dim> F_Inverse_T = ALETransformations 02911 ::get_F_Inverse_T<dim> (F_Inverse); 02912 02913 const double J = ALETransformations 02914 ::get_J<dim> (F); 02915 02916 const Tensor<2,dim> sigma_ALE = NSEALE 02917 ::get_stress_fluid_except_pressure_ALE<dim> 02918 (density_fluid, viscosity, grad_v, grad_v_T, F_Inverse, F_Inverse_T ); 02919 02920 Tensor<2,dim> stress_fluid; 02921 stress_fluid.clear(); 02922 stress_fluid = (J * sigma_ALE * F_Inverse_T); 02923 02924 Tensor<2,dim> fluid_pressure; 02925 fluid_pressure.clear(); 02926 fluid_pressure = (-pI * J * F_Inverse_T); 02927 02928 drag_lift_value -= (stress_fluid + fluid_pressure) * 02929 fe_face_values.normal_vector(q_point)* fe_face_values.JxW(q_point); 02930 } 02931 } 02932 } 02933 } 02934 02935 std::cout << "Drag: " << drag_lift_value[0] << std::endl; 02936 std::cout << "Lift: " << drag_lift_value[1] << std::endl; 02937 } 02938 02942 template<int dim> 02943 void FSIALEProblem<dim>::compute_functional_values() 02944 { 02945 double x1,y1; 02946 x1 = compute_point_value(Point<dim>(0.6,0.2), dim); 02947 y1 = compute_point_value(Point<dim>(0.6,0.2), dim+1); 02948 02949 std::cout << "------------------" << std::endl; 02950 std::cout << "DisX: " << x1 << std::endl; 02951 std::cout << "DisY: " << y1 << std::endl; 02952 std::cout << "------------------" << std::endl; 02953 02954 compute_drag_lift_fsi_fluid_tensor(); 02955 02956 std::cout << std::endl; 02957 } 02958 02968 template <int dim> 02969 void FSIALEProblem<dim>::run () 02970 { 02971 setup_system(); 02972 02973 std::cout << "\n==============================" 02974 << "=====================================" << std::endl; 02975 std::cout << "Parameters\n" 02976 << "==========\n" 02977 << "Density fluid: " << density_fluid << "\n" 02978 << "Density structure: " << density_structure << "\n" 02979 << "Viscosity fluid: " << viscosity << "\n" 02980 << "alpha_u: " << alpha_u << "\n" 02981 << "alpha_w: " << alpha_w << "\n" 02982 << "Lame coeff. mu: " << lame_coefficient_mu << "\n" 02983 << std::endl; 02984 02985 02986 const unsigned int output_skip = 5; 02987 do 02988 { 02989 std::cout << "Timestep " << timestep_number 02990 << " (" << time_stepping_scheme 02991 << ")" << ": " << time 02992 << " (" << timestep << ")" 02993 << "\n==============================" 02994 << "=====================================" 02995 << std::endl; 02996 02997 std::cout << std::endl; 02998 02999 // Compute next time step 03000 old_timestep_solution = solution; 03001 newton_iteration (time); 03002 time += timestep; 03003 03004 // Compute functional values: dx, dy, drag, lift 03005 std::cout << std::endl; 03006 compute_functional_values(); 03007 03008 // Write solutions 03009 if ((timestep_number % output_skip == 0)) 03010 output_results (timestep_number,solution); 03011 03012 03013 ++timestep_number; 03014 03015 } 03016 while (timestep_number <= max_no_timesteps); 03017 03018 03019 } 03020 03021 // The main function looks almost the same 03022 // as in all other deal.II tuturial steps. 03023 int main () 03024 { 03025 try 03026 { 03027 deallog.depth_console (0); 03028 03029 FSIALEProblem<2> fsi_problem(1); 03030 fsi_problem.run (); 03031 } 03032 catch (std::exception &exc) 03033 { 03034 std::cerr << std::endl << std::endl 03035 << "----------------------------------------------------" 03036 << std::endl; 03037 std::cerr << "Exception on processing: " << std::endl 03038 << exc.what() << std::endl 03039 << "Aborting!" << std::endl 03040 << "----------------------------------------------------" 03041 << std::endl; 03042 03043 return 1; 03044 } 03045 catch (...) 03046 { 03047 std::cerr << std::endl << std::endl 03048 << "----------------------------------------------------" 03049 << std::endl; 03050 std::cerr << "Unknown exception!" << std::endl 03051 << "Aborting!" << std::endl 03052 << "----------------------------------------------------" 03053 << std::endl; 03054 return 1; 03055 } 03056 03057 return 0; 03058 } 03059 03060 03061 03062
1.7.6.1