9 #ifndef __stokes_polynomial_h 10 #define __stokes_polynomial_h 12 #include <amandus/integrator.h> 13 #include <deal.II/base/polynomial.h> 14 #include <deal.II/base/tensor.h> 15 #include <deal.II/integrators/divergence.h> 16 #include <deal.II/integrators/l2.h> 17 #include <deal.II/integrators/laplace.h> 43 PolynomialRHS(
const Polynomials::Polynomial<double> curl_potential_1d,
44 const Polynomials::Polynomial<double> grad_potential_1d);
46 virtual void cell(DoFInfo<dim>& dinfo, IntegrationInfo<dim>& info)
const;
47 virtual void boundary(DoFInfo<dim>& dinfo, IntegrationInfo<dim>& info)
const;
48 virtual void face(DoFInfo<dim>& dinfo1, DoFInfo<dim>& dinfo2, IntegrationInfo<dim>& info1,
49 IntegrationInfo<dim>& info2)
const;
80 const Polynomials::Polynomial<double> grad_potential_1d);
82 virtual void cell(DoFInfo<dim>& dinfo, IntegrationInfo<dim>& info)
const;
83 virtual void boundary(DoFInfo<dim>& dinfo, IntegrationInfo<dim>& info)
const;
84 virtual void face(DoFInfo<dim>& dinfo1, DoFInfo<dim>& dinfo2, IntegrationInfo<dim>& info1,
85 IntegrationInfo<dim>& info2)
const;
129 PolynomialError(
const Polynomials::Polynomial<double> curl_potential_1d,
130 const Polynomials::Polynomial<double> grad_potential_1d);
132 virtual void cell(DoFInfo<dim>& dinfo, IntegrationInfo<dim>& info)
const;
133 virtual void boundary(DoFInfo<dim>& dinfo, IntegrationInfo<dim>& info)
const;
134 virtual void face(DoFInfo<dim>& dinfo1, DoFInfo<dim>& dinfo2, IntegrationInfo<dim>& info1,
135 IntegrationInfo<dim>& info2)
const;
146 const Polynomials::Polynomial<double> grad_potential_1d)
147 : curl_potential_1d(curl_potential_1d)
148 , grad_potential_1d(grad_potential_1d)
150 this->use_boundary =
false;
151 this->use_face =
false;
158 std::vector<std::vector<double>> rhs(dim,
159 std::vector<double>(info.fe_values(0).n_quadrature_points));
161 std::vector<double> px(4);
162 std::vector<double> py(4);
163 for (
unsigned int k = 0; k < info.fe_values(0).n_quadrature_points; ++k)
165 const double x = info.fe_values(0).quadrature_point(k)(0);
166 const double y = info.fe_values(0).quadrature_point(k)(1);
170 rhs[0][k] = -px[2] * py[1] - px[0] * py[3];
171 rhs[1][k] = px[3] * py[0] + px[1] * py[2];
178 rhs[0][k] += px[1] * py[0];
179 rhs[1][k] += px[0] * py[1];
182 L2::L2(dinfo.vector(0).block(0), info.fe_values(0), rhs);
194 IntegrationInfo<dim>&)
const 203 : curl_potential_1d(curl_potential_1d)
204 , grad_potential_1d(grad_potential_1d)
206 this->use_boundary =
true;
207 this->use_face =
true;
208 this->input_vector_names.push_back(
"Newton iterate");
215 Assert(info.values.size() >= 1, ExcDimensionMismatch(info.values.size(), 1));
216 Assert(info.gradients.size() >= 1, ExcDimensionMismatch(info.values.size(), 1));
218 std::vector<std::vector<double>> rhs(dim,
219 std::vector<double>(info.fe_values(0).n_quadrature_points));
221 std::vector<double> px(4);
222 std::vector<double> py(4);
223 for (
unsigned int k = 0; k < info.fe_values(0).n_quadrature_points; ++k)
225 const double x = info.fe_values(0).quadrature_point(k)(0);
226 const double y = info.fe_values(0).quadrature_point(k)(1);
230 rhs[0][k] = -px[2] * py[1] - px[0] * py[3];
231 rhs[1][k] = px[3] * py[0] + px[1] * py[2];
238 rhs[0][k] += px[1] * py[0];
239 rhs[1][k] += px[0] * py[1];
242 L2::L2(dinfo.vector(0).block(0), info.fe_values(0), rhs, -1.);
244 dinfo.vector(0).block(0), info.fe_values(0), make_slice(info.gradients[0], 0, dim));
245 Divergence::gradient_residual(
246 dinfo.vector(0).block(0), info.fe_values(0), info.values[0][dim], -1.);
248 dinfo.vector(0).block(1), info.fe_values(1), make_slice(info.gradients[0], 0, dim), 1.);
255 std::vector<std::vector<double>> null(
256 dim, std::vector<double>(info.fe_values(0).n_quadrature_points, 0.));
258 const unsigned int deg = info.fe_values(0).get_fe().tensor_degree();
259 Laplace::nitsche_residual(dinfo.vector(0).block(0),
261 make_slice(info.values[0], 0, dim),
262 make_slice(info.gradients[0], 0, dim),
264 Laplace::compute_penalty(dinfo, dinfo, deg, deg));
270 IntegrationInfo<dim>& info1, IntegrationInfo<dim>& info2)
const 272 const unsigned int deg = info1.fe_values(0).get_fe().tensor_degree();
273 Laplace::ip_residual(dinfo1.vector(0).block(0),
274 dinfo2.vector(0).block(0),
277 make_slice(info1.values[0], 0, dim),
278 make_slice(info1.gradients[0], 0, dim),
279 make_slice(info2.values[0], 0, dim),
280 make_slice(info2.gradients[0], 0, dim),
281 Laplace::compute_penalty(dinfo1, dinfo2, deg, deg));
289 : curl_potential_1d(curl_potential_1d)
290 , grad_potential_1d(grad_potential_1d)
306 this->use_boundary =
false;
307 this->use_face =
false;
316 std::vector<double> px(3);
317 std::vector<double> py(3);
318 for (
unsigned int k = 0; k < info.fe_values(0).n_quadrature_points; ++k)
320 const double x = info.fe_values(0).quadrature_point(k)(0);
321 const double y = info.fe_values(0).quadrature_point(k)(1);
324 const double dx = info.fe_values(0).JxW(k);
326 Tensor<1, dim> Du0 = info.gradients[0][0][k];
327 Du0[0] -= px[1] * py[1];
328 Du0[1] -= px[0] * py[2];
329 Tensor<1, dim> Du1 = info.gradients[0][1][k];
330 Du1[0] += px[2] * py[0];
331 Du1[1] += px[1] * py[1];
332 double u0 = info.values[0][0][k];
334 double u1 = info.values[0][1][k];
339 double p = info.values[0][dim][k];
341 Tensor<1, dim> Dp = info.gradients[0][dim][k];
342 Dp[0] += px[1] * py[0];
343 Dp[1] += px[0] * py[1];
345 dinfo.value(0) += (u0 * u0 + u1 * u1) * dx;
347 dinfo.value(1) += ((Du0 * Du0) + (Du1 * Du1)) * dx;
349 dinfo.value(2) = Divergence::norm(info.fe_values(0), make_slice(info.gradients[0], 0, dim));
351 dinfo.value(3) += p * p * dx;
353 dinfo.value(4) = std::max(dinfo.value(4), std::abs(p));
355 dinfo.value(5) += Dp * Dp * dx;
357 dinfo.value(6) += p * dx;
360 for (
unsigned int i = 0; i <= 4; ++i)
362 dinfo.value(i) = std::sqrt(dinfo.value(i));
374 IntegrationInfo<dim>&)
const virtual void boundary(DoFInfo< dim > &dinfo, IntegrationInfo< dim > &info) const
Definition: stokes/polynomial.h:187
virtual void boundary(DoFInfo< dim > &dinfo, IntegrationInfo< dim > &info) const
Definition: stokes/polynomial.h:367
std::vector< std::string > error_names
Definition: integrator.h:120
PolynomialResidual(const Polynomials::Polynomial< double > curl_potential_1d, const Polynomials::Polynomial< double > grad_potential_1d)
Definition: stokes/polynomial.h:201
Polynomials::Polynomial< double > grad_potential_1d
Definition: stokes/polynomial.h:89
Polynomials::Polynomial< double > grad_potential_1d
Definition: stokes/polynomial.h:139
Polynomials::Polynomial< double > curl_potential_1d
Definition: stokes/polynomial.h:88
Definition: stokes/polynomial.h:126
virtual void boundary(DoFInfo< dim > &dinfo, IntegrationInfo< dim > &info) const
Definition: stokes/polynomial.h:253
virtual void face(DoFInfo< dim > &dinfo1, DoFInfo< dim > &dinfo2, IntegrationInfo< dim > &info1, IntegrationInfo< dim > &info2) const
Definition: stokes/polynomial.h:193
Definition: stokes/polynomial.h:40
virtual void cell(DoFInfo< dim > &dinfo, IntegrationInfo< dim > &info) const
Definition: stokes/polynomial.h:213
Definition: stokes/eigen.h:37
virtual void cell(DoFInfo< dim > &dinfo, IntegrationInfo< dim > &info) const
Definition: stokes/polynomial.h:312
virtual void face(DoFInfo< dim > &dinfo1, DoFInfo< dim > &dinfo2, IntegrationInfo< dim > &info1, IntegrationInfo< dim > &info2) const
Definition: stokes/polynomial.h:269
virtual void cell(DoFInfo< dim > &dinfo, IntegrationInfo< dim > &info) const
Definition: stokes/polynomial.h:156
Polynomials::Polynomial< double > grad_potential_1d
Definition: stokes/polynomial.h:53
std::vector< unsigned int > error_types
Definition: integrator.h:113
PolynomialError(const Polynomials::Polynomial< double > curl_potential_1d, const Polynomials::Polynomial< double > grad_potential_1d)
Definition: stokes/polynomial.h:287
void cell_residual(dealii::Vector< number > &result, const dealii::FEValuesBase< dim > &fe, const dealii::VectorSlice< const std::vector< std::vector< dealii::Tensor< 1, dim >>>> &input, double lambda=0., double mu=1.)
Definition: elasticity/integrators.h:23
Polynomials::Polynomial< double > curl_potential_1d
Definition: stokes/polynomial.h:138
Definition: stokes/polynomial.h:76
Polynomials::Polynomial< double > curl_potential_1d
Definition: stokes/polynomial.h:52
virtual void face(DoFInfo< dim > &dinfo1, DoFInfo< dim > &dinfo2, IntegrationInfo< dim > &info1, IntegrationInfo< dim > &info2) const
Definition: stokes/polynomial.h:373
unsigned int error_type(unsigned int i) const
Definition: integrator.h:437
Definition: integrator.h:29