//---------------------------------------------------------------------------
// $Id$
//
-// Copyright (C) 2010, 2011, 2012 by the deal.II authors
+// Copyright (C) 2010, 2011, 2012, 2013 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
namespace Elasticity
{
/**
- * Scalar product of symmetric gradients.
+ * The linear elasticity operator in weak form, namely double
+ * contraction of symmetric gradients.
*
* \f[
- * (\varepsilon(u), \varepsilon(v))
+ * \int_Z \varepsilon(u): \varepsilon(v)\,dx
* \f]
*/
template <int dim>
/**
- * The weak boundary condition
- * of Nitsche type for
- * symmetric gradients.
+ * Vector-valued residual operator for linear elasticity in weak form
+ *
+ * \f[
+ * - \int_Z \varepsilon(u): \varepsilon(v) \,dx
+ * \f]
+ */
+ template <int dim>
+ inline void
+ cell_residual (
+ Vector<double> &result,
+ const FEValuesBase<dim> &fe,
+ const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &input,
+ double factor = 1.)
+ {
+ const unsigned int nq = fe.n_quadrature_points;
+ const unsigned int n_dofs = fe.dofs_per_cell;
+ AssertDimension(fe.get_fe().n_components(), dim);
+
+ AssertVectorVectorDimension(input, dim, fe.n_quadrature_points);
+ Assert(result.size() == n_dofs, ExcDimensionMismatch(result.size(), n_dofs));
+
+ for (unsigned int k=0; k<nq; ++k)
+ {
+ const double dx = factor * fe.JxW(k);
+ for (unsigned int i=0; i<n_dofs; ++i)
+ for (unsigned int d1=0; d1<dim; ++d1)
+ for (unsigned int d2=0; d2<dim; ++d2)
+ {
+ result(i) += dx * .25 *
+ (input[d1][k][d2] + input[d2][k][d1]) *
+ (fe.shape_grad_component(i,k,d1)[d2] + fe.shape_grad_component(i,k,d2)[d1]);
+ }
+ }
+ }
+
+
+ /**
+ * The weak boundary condition of Nitsche type for linear
+ * elasticity:
+ * @f[
+ * \int_F \Bigl(\gamma (u-g) \cdot v - n^T \epsilon(u) v - (u-g) \epsilon(v) n^T\Bigr)\;ds.
+ * @f]
*/
template <int dim>
inline void nitsche_matrix (
}
}
+ /**
+ * Weak boundary condition for the elasticity operator by Nitsche,
+ * namely on the face <i>F</i> the vector
+ * @f[
+ * \int_F \Bigl(\gamma (u-g) \cdot v - n^T \epsilon(u) v - (u-g) \epsilon(v) n^T\Bigr)\;ds.
+ * @f]
+ *
+ * Here, <i>u</i> is the finite element function whose values and
+ * gradient are given in the arguments <tt>input</tt> and
+ * <tt>Dinput</tt>, respectively. <i>g</i> is the inhomogeneous
+ * boundary value in the argument <tt>data</tt>. $n$ is the outer
+ * normal vector and $\gamma$ is the usual penalty parameter.
+ *
+ * @author Guido Kanschat
+ * @date 2013
+ */
+ template <int dim>
+ void nitsche_residual (
+ Vector<double> &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.)
+ {
+ const unsigned int n_dofs = fe.dofs_per_cell;
+ AssertVectorVectorDimension(input, dim, fe.n_quadrature_points);
+ AssertVectorVectorDimension(Dinput, dim, fe.n_quadrature_points);
+ AssertVectorVectorDimension(data, dim, fe.n_quadrature_points);
+
+ for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
+ {
+ const double dx = factor * fe.JxW(k);
+ const Point<dim> &n = fe.normal_vector(k);
+ for (unsigned int i=0; i<n_dofs; ++i)
+ for (unsigned int d1=0; d1<dim; ++d1)
+ {
+ const double u= input[d1][k];
+ const double v= fe.shape_value_component(i,k,d1);
+ const double g= data[d1][k];
+ result(i) += dx + 2.*penalty * (u-g) * v;
+
+ for (unsigned int d2=0; d2<dim; ++d2)
+ {
+ // v . nabla u n
+ result(i) -= .5*dx* v * Dinput[d1][k][d2] * n(d2);
+ // v . (nabla u)^T n
+ result(i) -= .5*dx* v * Dinput[d2][k][d1] * n(d2);
+ // u nabla v n
+ result(i) -= .5*dx * (u-g) * fe.shape_grad_component(i,k,d1)[d2] * n(d2);
+ // u (nabla v)^T n
+ result(i) -= .5*dx * (u-g) * fe.shape_grad_component(i,k,d2)[d1] * n(d2);
+ }
+ }
+ }
+ }
+
/**
* The interior penalty flux
* for symmetric gradients.
}
}
}
+ /**
+ * Elasticity residual term for the symmetric interior penalty method.
+ *
+ * @author Guido Kanschat
+ * @date 2013
+ */
+ template<int dim>
+ void
+ ip_residual(
+ Vector<double> &result1,
+ Vector<double> &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.)
+ {
+ const unsigned int n1 = fe1.dofs_per_cell;
+
+ AssertDimension(fe1.get_fe().n_components(), dim);
+ AssertDimension(fe2.get_fe().n_components(), dim);
+ AssertVectorVectorDimension(input1, dim, fe1.n_quadrature_points);
+ AssertVectorVectorDimension(Dinput1, dim, fe1.n_quadrature_points);
+ AssertVectorVectorDimension(input2, dim, fe2.n_quadrature_points);
+ AssertVectorVectorDimension(Dinput2, dim, fe2.n_quadrature_points);
+
+ const double nu1 = int_factor;
+ const double nu2 = (ext_factor < 0) ? int_factor : ext_factor;
+ const double penalty = .5 * pen * (nu1 + nu2);
+
+
+ for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
+ {
+ const double dx = fe1.JxW(k);
+ const Point<dim> &n = fe1.normal_vector(k);
+
+ for (unsigned int i=0; i<n1; ++i)
+ for (unsigned int d1=0; d1<dim; ++d1)
+ {
+ const double v1 = fe1.shape_value_component(i,k,d1);
+ const double v2 = fe2.shape_value_component(i,k,d1);
+ const double u1 = input1[d1][k];
+ const double u2 = input2[d1][k];
+
+ result1(i) += dx * penalty * u1*v1;
+ result1(i) -= dx * penalty * u2*v1;
+ result2(i) -= dx * penalty * u1*v2;
+ result2(i) += dx * penalty * u2*v2;
+
+ for (unsigned int d2=0; d2<dim; ++d2)
+ {
+ // v . nabla u n
+ result1(i) -= .25*dx* (nu1*Dinput1[d1][k][d2]+nu2*Dinput1[d1][k][d2]) * n(d2) * v1;
+ result2(i) += .25*dx* (nu1*Dinput1[d1][k][d2]+nu2*Dinput1[d1][k][d2]) * n(d2) * v1;
+ // v . (nabla u)^T n
+ result1(i) -= .25*dx* (nu1*Dinput1[d2][k][d1]+nu2*Dinput1[d2][k][d1]) * n(d2) * v1;
+ result2(i) += .25*dx* (nu1*Dinput1[d2][k][d1]+nu2*Dinput1[d2][k][d1]) * n(d2) * v1;
+ // u nabla v n
+ result1(i) -= .25*dx* nu1*fe1.shape_grad_component(i,k,d1)[d2] * n(d2) * (u1-u2);
+ result2(i) -= .25*dx* nu2*fe2.shape_grad_component(i,k,d1)[d2] * n(d2) * (u1-u2);
+ // u (nabla v)^T n
+ result1(i) -= .25*dx* nu1*fe1.shape_grad_component(i,k,d2)[d1] * n(d2) * (u1-u2);
+ result2(i) -= .25*dx* nu2*fe2.shape_grad_component(i,k,d2)[d1] * n(d2) * (u1-u2);
+ }
+ }
+ }
+ }
}
}
projects / dealii-svn.git / commitdiff