00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__integrators_elasticity_h
00013 #define __deal2__integrators_elasticity_h
00014
00015
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/exceptions.h>
00018 #include <deal.II/base/quadrature.h>
00019 #include <deal.II/lac/full_matrix.h>
00020 #include <deal.II/fe/mapping.h>
00021 #include <deal.II/fe/fe_values.h>
00022 #include <deal.II/meshworker/dof_info.h>
00023
00024 DEAL_II_NAMESPACE_OPEN
00025
00026 namespace LocalIntegrators
00027 {
00035 namespace Elasticity
00036 {
00044 template <int dim>
00045 inline void cell_matrix (
00046 FullMatrix<double>& M,
00047 const FEValuesBase<dim>& fe,
00048 const double factor = 1.)
00049 {
00050 const unsigned int n_dofs = fe.dofs_per_cell;
00051
00052 AssertDimension(fe.get_fe().n_components(), dim);
00053 AssertDimension(M.m(), n_dofs);
00054 AssertDimension(M.n(), n_dofs);
00055
00056 for (unsigned k=0;k<fe.n_quadrature_points;++k)
00057 {
00058 const double dx = factor * fe.JxW(k);
00059 for (unsigned i=0;i<n_dofs;++i)
00060 for (unsigned j=0;j<n_dofs;++j)
00061 for (unsigned int d1=0;d1<dim;++d1)
00062 for (unsigned int d2=0;d2<dim;++d2)
00063 M(i,j) += dx * .25 *
00064 (fe.shape_grad_component(j,k,d1)[d2] + fe.shape_grad_component(j,k,d2)[d1]) *
00065 (fe.shape_grad_component(i,k,d1)[d2] + fe.shape_grad_component(i,k,d2)[d1]);
00066 }
00067 }
00068
00069
00075 template <int dim>
00076 inline void nitsche_matrix (
00077 FullMatrix<double>& M,
00078 const FEValuesBase<dim>& fe,
00079 double penalty,
00080 double factor = 1.)
00081 {
00082 const unsigned int n_dofs = fe.dofs_per_cell;
00083
00084 AssertDimension(fe.get_fe().n_components(), dim);
00085 AssertDimension(M.m(), n_dofs);
00086 AssertDimension(M.n(), n_dofs);
00087
00088 for (unsigned k=0;k<fe.n_quadrature_points;++k)
00089 {
00090 const double dx = factor * fe.JxW(k);
00091 const Point<dim>& n = fe.normal_vector(k);
00092 for (unsigned i=0;i<n_dofs;++i)
00093 for (unsigned j=0;j<n_dofs;++j)
00094 for (unsigned int d1=0;d1<dim;++d1)
00095 {
00096 const double u = fe.shape_value_component(j,k,d1);
00097 const double v = fe.shape_value_component(i,k,d1);
00098 M(i,j) += dx * penalty * u * v;
00099 for (unsigned int d2=0;d2<dim;++d2)
00100 {
00101
00102 M(i,j) -= .5*dx* fe.shape_grad_component(j,k,d1)[d2] *n(d2)* v;
00103
00104 M(i,j) -= .5*dx* fe.shape_grad_component(j,k,d2)[d1] *n(d2)* v;
00105
00106 M(i,j) -= .5*dx* fe.shape_grad_component(i,k,d1)[d2] *n(d2)* u;
00107
00108 M(i,j) -= .5*dx* fe.shape_grad_component(i,k,d2)[d1] *n(d2)* u;
00109 }
00110 }
00111 }
00112 }
00113
00118 template <int dim>
00119 inline void ip_matrix (
00120 FullMatrix<double>& M11,
00121 FullMatrix<double>& M12,
00122 FullMatrix<double>& M21,
00123 FullMatrix<double>& M22,
00124 const FEValuesBase<dim>& fe1,
00125 const FEValuesBase<dim>& fe2,
00126 const double pen,
00127 const double int_factor = 1.,
00128 const double ext_factor = -1.)
00129 {
00130 const unsigned int n_dofs = fe1.dofs_per_cell;
00131
00132 AssertDimension(fe1.get_fe().n_components(), dim);
00133 AssertDimension(fe2.get_fe().n_components(), dim);
00134 AssertDimension(M11.m(), n_dofs);
00135 AssertDimension(M11.n(), n_dofs);
00136 AssertDimension(M12.m(), n_dofs);
00137 AssertDimension(M12.n(), n_dofs);
00138 AssertDimension(M21.m(), n_dofs);
00139 AssertDimension(M21.n(), n_dofs);
00140 AssertDimension(M22.m(), n_dofs);
00141 AssertDimension(M22.n(), n_dofs);
00142
00143 const double nu1 = int_factor;
00144 const double nu2 = (ext_factor < 0) ? int_factor : ext_factor;
00145 const double penalty = .5 * pen * (nu1 + nu2);
00146
00147 for (unsigned k=0;k<fe1.n_quadrature_points;++k)
00148 {
00149 const double dx = fe1.JxW(k);
00150 const Point<dim>& n = fe1.normal_vector(k);
00151 for (unsigned i=0;i<n_dofs;++i)
00152 for (unsigned j=0;j<n_dofs;++j)
00153 for (unsigned int d1=0;d1<dim;++d1)
00154 {
00155 const double u1 = fe1.shape_value_component(j,k,d1);
00156 const double u2 = fe2.shape_value_component(j,k,d1);
00157 const double v1 = fe1.shape_value_component(i,k,d1);
00158 const double v2 = fe2.shape_value_component(i,k,d1);
00159
00160 M11(i,j) += dx * penalty * u1*v1;
00161 M12(i,j) -= dx * penalty * u2*v1;
00162 M21(i,j) -= dx * penalty * u1*v2;
00163 M22(i,j) += dx * penalty * u2*v2;
00164
00165 for (unsigned int d2=0;d2<dim;++d2)
00166 {
00167
00168 M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(j,k,d1)[d2] * n(d2) * v1;
00169 M12(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(j,k,d1)[d2] * n(d2) * v1;
00170 M21(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(j,k,d1)[d2] * n(d2) * v2;
00171 M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(j,k,d1)[d2] * n(d2) * v2;
00172
00173 M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(j,k,d2)[d1] * n(d2) * v1;
00174 M12(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(j,k,d2)[d1] * n(d2) * v1;
00175 M21(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(j,k,d2)[d1] * n(d2) * v2;
00176 M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(j,k,d2)[d1] * n(d2) * v2;
00177
00178 M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(i,k,d1)[d2] * n(d2) * u1;
00179 M12(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(i,k,d1)[d2] * n(d2) * u2;
00180 M21(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(i,k,d1)[d2] * n(d2) * u1;
00181 M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(i,k,d1)[d2] * n(d2) * u2;
00182
00183 M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(i,k,d2)[d1] * n(d2) * u1;
00184 M12(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(i,k,d2)[d1] * n(d2) * u2;
00185 M21(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(i,k,d2)[d1] * n(d2) * u1;
00186 M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(i,k,d2)[d1] * n(d2) * u2;
00187 }
00188 }
00189 }
00190 }
00191 }
00192 }
00193
00194 DEAL_II_NAMESPACE_CLOSE
00195
00196 #endif
00197