17 #ifndef __deal2__integrators_elasticity_h
18 #define __deal2__integrators_elasticity_h
21 #include <deal.II/base/config.h>
23 #include <deal.II/base/quadrature.h>
24 #include <deal.II/lac/full_matrix.h>
25 #include <deal.II/fe/mapping.h>
26 #include <deal.II/fe/fe_values.h>
27 #include <deal.II/meshworker/dof_info.h>
54 const double factor = 1.)
64 const double dx = factor * fe.
JxW(k);
65 for (
unsigned int i=0; i<n_dofs; ++i)
66 for (
unsigned int j=0; j<n_dofs; ++j)
67 for (
unsigned int d1=0; d1<dim; ++d1)
68 for (
unsigned int d2=0; d2<dim; ++d2)
83 template <
int dim,
typename number>
98 for (
unsigned int k=0; k<nq; ++k)
100 const double dx = factor * fe.
JxW(k);
101 for (
unsigned int i=0; i<n_dofs; ++i)
102 for (
unsigned int d1=0; d1<dim; ++d1)
103 for (
unsigned int d2=0; d2<dim; ++d2)
105 result(i) += dx * .25 *
106 (input[d1][k][d2] + input[d2][k][d1]) *
135 const double dx = factor * fe.
JxW(k);
137 for (
unsigned int i=0; i<n_dofs; ++i)
138 for (
unsigned int j=0; j<n_dofs; ++j)
139 for (
unsigned int d1=0; d1<dim; ++d1)
143 M(i,j) += dx * 2. * penalty * u * v;
144 for (
unsigned int d2=0; d2<dim; ++d2)
175 template <
int dim,
typename number>
179 const VectorSlice<
const std::vector<std::vector<double> > > &input,
181 const VectorSlice<
const std::vector<std::vector<double> > > &data,
192 const double dx = factor * fe.
JxW(k);
194 for (
unsigned int i=0; i<n_dofs; ++i)
195 for (
unsigned int d1=0; d1<dim; ++d1)
197 const double u= input[d1][k];
199 const double g= data[d1][k];
200 result(i) += dx * 2.*penalty * (u-g) * v;
202 for (
unsigned int d2=0; d2<dim; ++d2)
205 result(i) -= .5*dx* v * Dinput[d1][k][d2] * n(d2);
207 result(i) -= .5*dx* v * Dinput[d2][k][d1] * n(d2);
232 template <
int dim,
typename number>
236 const VectorSlice<
const std::vector<std::vector<double> > > &input,
247 const double dx = factor * fe.
JxW(k);
249 for (
unsigned int i=0; i<n_dofs; ++i)
250 for (
unsigned int d1=0; d1<dim; ++d1)
252 const double u= input[d1][k];
254 result(i) += dx * 2.*penalty * u * v;
256 for (
unsigned int d2=0; d2<dim; ++d2)
259 result(i) -= .5*dx* v * Dinput[d1][k][d2] * n(d2);
261 result(i) -= .5*dx* v * Dinput[d2][k][d1] * n(d2);
284 const double int_factor = 1.,
285 const double ext_factor = -1.)
300 const double nu1 = int_factor;
301 const double nu2 = (ext_factor < 0) ? int_factor : ext_factor;
302 const double penalty = .5 * pen * (nu1 + nu2);
306 const double dx = fe1.
JxW(k);
308 for (
unsigned int i=0; i<n_dofs; ++i)
309 for (
unsigned int j=0; j<n_dofs; ++j)
310 for (
unsigned int d1=0; d1<dim; ++d1)
317 M11(i,j) += dx * penalty * u1*v1;
318 M12(i,j) -= dx * penalty * u2*v1;
319 M21(i,j) -= dx * penalty * u1*v2;
320 M22(i,j) += dx * penalty * u2*v2;
322 for (
unsigned int d2=0; d2<dim; ++d2)
354 template<
int dim,
typename number>
361 const VectorSlice<
const std::vector<std::vector<double> > > &input1,
363 const VectorSlice<
const std::vector<std::vector<double> > > &input2,
366 double int_factor = 1.,
367 double ext_factor = -1.)
378 const double nu1 = int_factor;
379 const double nu2 = (ext_factor < 0) ? int_factor : ext_factor;
380 const double penalty = .5 * pen * (nu1 + nu2);
385 const double dx = fe1.
JxW(k);
388 for (
unsigned int i=0; i<n1; ++i)
389 for (
unsigned int d1=0; d1<dim; ++d1)
393 const double u1 = input1[d1][k];
394 const double u2 = input2[d1][k];
396 result1(i) += dx * penalty * u1*v1;
397 result1(i) -= dx * penalty * u2*v1;
398 result2(i) -= dx * penalty * u1*v2;
399 result2(i) += dx * penalty * u2*v2;
401 for (
unsigned int d2=0; d2<dim; ++d2)
404 result1(i) -= .25*dx* (nu1*Dinput1[d1][k][d2]+nu2*Dinput2[d1][k][d2]) * n(d2) * v1;
405 result2(i) += .25*dx* (nu1*Dinput1[d1][k][d2]+nu2*Dinput2[d1][k][d2]) * n(d2) * v2;
407 result1(i) -= .25*dx* (nu1*Dinput1[d2][k][d1]+nu2*Dinput2[d2][k][d1]) * n(d2) * v1;
408 result2(i) += .25*dx* (nu1*Dinput1[d2][k][d1]+nu2*Dinput2[d2][k][d1]) * n(d2) * v2;
422 DEAL_II_NAMESPACE_CLOSE
double shape_value_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
void nitsche_residual_homogeneous(Vector< number > &result, const FEValuesBase< dim > &fe, const VectorSlice< const std::vector< std::vector< double > > > &input, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim > > > > &Dinput, double penalty, double factor=1.)
#define AssertDimension(dim1, dim2)
Tensor< 1, spacedim > shape_grad_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
const unsigned int dofs_per_cell
const Point< spacedim > & normal_vector(const unsigned int i) const
#define AssertVectorVectorDimension(vec, dim1, dim2)
Library of integrals over cells and faces.
void cell_residual(Vector< number > &result, const FEValuesBase< dim > &fe, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim > > > > &input, double factor=1.)
#define Assert(cond, exc)
void ip_matrix(FullMatrix< double > &M11, FullMatrix< double > &M12, FullMatrix< double > &M21, FullMatrix< double > &M22, const FEValuesBase< dim > &fe1, const FEValuesBase< dim > &fe2, const double pen, const double int_factor=1., const double ext_factor=-1.)
const unsigned int n_quadrature_points
void nitsche_residual(Vector< number > &result, const FEValuesBase< dim > &fe, const VectorSlice< const std::vector< std::vector< double > > > &input, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim > > > > &Dinput, const VectorSlice< const std::vector< std::vector< double > > > &data, double penalty, double factor=1.)
void nitsche_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
const FiniteElement< dim, spacedim > & get_fe() const
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
double JxW(const unsigned int quadrature_point) const
void ip_residual(Vector< number > &result1, Vector< number > &result2, const FEValuesBase< dim > &fe1, const FEValuesBase< dim > &fe2, const VectorSlice< const std::vector< std::vector< double > > > &input1, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim > > > > &Dinput1, const VectorSlice< const std::vector< std::vector< double > > > &input2, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim > > > > &Dinput2, double pen, double int_factor=1., double ext_factor=-1.)
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)