00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__integrators_laplace_h
00013 #define __deal2__integrators_laplace_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 Laplace
00036 {
00051 template<int dim>
00052 void cell_matrix (
00053 FullMatrix<double>& M,
00054 const FEValuesBase<dim>& fe,
00055 const double factor = 1.)
00056 {
00057 const unsigned int n_dofs = fe.dofs_per_cell;
00058 const unsigned int n_components = fe.get_fe().n_components();
00059
00060 for (unsigned k=0;k<fe.n_quadrature_points;++k)
00061 {
00062 const double dx = fe.JxW(k) * factor;
00063 for (unsigned i=0;i<n_dofs;++i)
00064 {
00065 for (unsigned j=0;j<n_dofs;++j)
00066 for (unsigned int d=0;d<n_components;++d)
00067 M(i,j) += dx *
00068 (fe.shape_grad_component(j,k,d) * fe.shape_grad_component(i,k,d));
00069 }
00070 }
00071 }
00085 template <int dim>
00086 void nitsche_matrix (
00087 FullMatrix<double>& M,
00088 const FEValuesBase<dim>& fe,
00089 double penalty,
00090 double factor = 1.)
00091 {
00092 const unsigned int n_dofs = fe.dofs_per_cell;
00093 const unsigned int n_comp = fe.get_fe().n_components();
00094
00095 Assert (M.m() == n_dofs, ExcDimensionMismatch(M.m(), n_dofs));
00096 Assert (M.n() == n_dofs, ExcDimensionMismatch(M.n(), n_dofs));
00097
00098 for (unsigned k=0;k<fe.n_quadrature_points;++k)
00099 {
00100 const double dx = fe.JxW(k) * factor;
00101 const Point<dim>& n = fe.normal_vector(k);
00102 for (unsigned i=0;i<n_dofs;++i)
00103 for (unsigned j=0;j<n_dofs;++j)
00104 for (unsigned int d=0;d<n_comp;++d)
00105 M(i,j) += dx *
00106 (2. * fe.shape_value_component(i,k,d) * penalty * fe.shape_value_component(j,k,d)
00107 - (n * fe.shape_grad_component(i,k,d)) * fe.shape_value_component(j,k,d)
00108 - (n * fe.shape_grad_component(j,k,d)) * fe.shape_value_component(i,k,d));
00109 }
00110 }
00111
00130 template <int dim>
00131 void nitsche_residual (
00132 Vector<double>& result,
00133 const FEValuesBase<dim>& fe,
00134 const VectorSlice<const std::vector<std::vector<double> > >& input,
00135 const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > >& Dinput,
00136 const VectorSlice<const std::vector<std::vector<double> > >& data,
00137 double penalty,
00138 double factor = 1.)
00139 {
00140 const unsigned int n_dofs = fe.dofs_per_cell;
00141
00142 const unsigned int n_comp = fe.get_fe().n_components();
00143 AssertVectorVectorDimension(input, n_comp, fe.n_quadrature_points);
00144 AssertVectorVectorDimension(Dinput, n_comp, fe.n_quadrature_points);
00145 AssertVectorVectorDimension(data, n_comp, fe.n_quadrature_points);
00146
00147 for (unsigned k=0;k<fe.n_quadrature_points;++k)
00148 {
00149 const double dx = factor * fe.JxW(k);
00150 const Point<dim>& n = fe.normal_vector(k);
00151 for (unsigned i=0;i<n_dofs;++i)
00152 for (unsigned int d=0;d<n_comp;++d)
00153 {
00154 const double dnv = fe.shape_grad_component(i,k,d) * n;
00155 const double dnu = Dinput[d][k] * n;
00156 const double v= fe.shape_value_component(i,k,d);
00157 const double u= input[d][k];
00158 const double g= data[d][k];
00159
00160 result(i) += dx*(2.*penalty*(u-g)*v - dnv*(u-g) - dnu*v);
00161 }
00162 }
00163 }
00164
00178 template <int dim>
00179 void nitsche_residual (
00180 Vector<double>& result,
00181 const FEValuesBase<dim>& fe,
00182 const std::vector<double>& input,
00183 const std::vector<Tensor<1,dim> >& Dinput,
00184 const std::vector<double>& data,
00185 double penalty,
00186 double factor = 1.)
00187 {
00188 const unsigned int n_dofs = fe.dofs_per_cell;
00189 AssertDimension(input.size(), fe.n_quadrature_points);
00190 AssertDimension(Dinput.size(), fe.n_quadrature_points);
00191 AssertDimension(data.size(), fe.n_quadrature_points);
00192
00193 for (unsigned k=0;k<fe.n_quadrature_points;++k)
00194 {
00195 const double dx = factor * fe.JxW(k);
00196 const Point<dim>& n = fe.normal_vector(k);
00197 for (unsigned i=0;i<n_dofs;++i)
00198 {
00199 const double dnv = fe.shape_grad(i,k) * n;
00200 const double dnu = Dinput[k] * n;
00201 const double v= fe.shape_value(i,k);
00202 const double u= input[k];
00203 const double g= data[k];
00204
00205 result(i) += dx*(2.*penalty*(u-g)*v - dnv*(u-g) - dnu*v);
00206 }
00207 }
00208 }
00209
00230 template <int dim>
00231 void ip_matrix (
00232 FullMatrix<double>& M11,
00233 FullMatrix<double>& M12,
00234 FullMatrix<double>& M21,
00235 FullMatrix<double>& M22,
00236 const FEValuesBase<dim>& fe1,
00237 const FEValuesBase<dim>& fe2,
00238 double penalty,
00239 double factor1 = 1.,
00240 double factor2 = -1.)
00241 {
00242 const unsigned int n_dofs = fe1.dofs_per_cell;
00243 AssertDimension(M11.n(), n_dofs);
00244 AssertDimension(M11.m(), n_dofs);
00245 AssertDimension(M12.n(), n_dofs);
00246 AssertDimension(M12.m(), n_dofs);
00247 AssertDimension(M21.n(), n_dofs);
00248 AssertDimension(M21.m(), n_dofs);
00249 AssertDimension(M22.n(), n_dofs);
00250 AssertDimension(M22.m(), n_dofs);
00251
00252 const double nui = factor1;
00253 const double nue = (factor2 < 0) ? factor1 : factor2;
00254
00255 for (unsigned k=0;k<fe1.n_quadrature_points;++k)
00256 {
00257 const double dx = fe1.JxW(k);
00258 const Point<dim>& n = fe1.normal_vector(k);
00259 for (unsigned int d=0;d<fe1.get_fe().n_components();++d)
00260 {
00261 for (unsigned i=0;i<n_dofs;++i)
00262 {
00263 for (unsigned j=0;j<n_dofs;++j)
00264 {
00265 const double vi = fe1.shape_value_component(i,k,d);
00266 const double dnvi = n * fe1.shape_grad_component(i,k,d);
00267 const double ve = fe2.shape_value_component(i,k,d);
00268 const double dnve = n * fe2.shape_grad_component(i,k,d);
00269 const double ui = fe1.shape_value_component(j,k,d);
00270 const double dnui = n * fe1.shape_grad_component(j,k,d);
00271 const double ue = fe2.shape_value_component(j,k,d);
00272 const double dnue = n * fe2.shape_grad_component(j,k,d);
00273 M11(i,j) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+penalty*ui*vi);
00274 M12(i,j) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-penalty*vi*ue);
00275 M21(i,j) += dx*(-.5*nue*dnve*ui+.5*nui*dnui*ve-penalty*ui*ve);
00276 M22(i,j) += dx*( .5*nue*dnve*ue+.5*nue*dnue*ve+penalty*ue*ve);
00277 }
00278 }
00279 }
00280 }
00281 }
00282
00297 template <int dim>
00298 double compute_penalty(
00299 const MeshWorker::DoFInfo<dim>& dinfo1,
00300 const MeshWorker::DoFInfo<dim>& dinfo2,
00301 unsigned int deg1,
00302 unsigned int deg2)
00303 {
00304 const unsigned int normal1 = GeometryInfo<dim>::unit_normal_direction[dinfo1.face_number];
00305 const unsigned int normal2 = GeometryInfo<dim>::unit_normal_direction[dinfo2.face_number];
00306 const unsigned int deg1sq = (deg1 == 0) ? 1 : deg1 * (deg1+1);
00307 const unsigned int deg2sq = (deg2 == 0) ? 1 : deg2 * (deg2+1);
00308
00309 double penalty1 = deg1sq / dinfo1.cell->extent_in_direction(normal1);
00310 double penalty2 = deg2sq / dinfo2.cell->extent_in_direction(normal2);
00311 if (dinfo1.cell->has_children() ^ dinfo2.cell->has_children())
00312 {
00313 Assert (dinfo1.face == dinfo2.face, ExcInternalError());
00314 Assert (dinfo1.face->has_children(), ExcInternalError());
00315 penalty1 *= 2;
00316 }
00317 const double penalty = 0.5*(penalty1 + penalty2);
00318 return penalty;
00319 }
00320 }
00321 }
00322
00323
00324 DEAL_II_NAMESPACE_CLOSE
00325
00326 #endif
00327