From: Wolfgang Bangerth Date: Fri, 2 Mar 2018 17:12:44 +0000 (-0700) Subject: Merge the two places where we compute beta in step-12 into one function. X-Git-Tag: v9.0.0-rc1~366^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F5988%2Fhead;p=dealii.git Merge the two places where we compute beta in step-12 into one function. --- diff --git a/examples/step-12/step-12.cc b/examples/step-12/step-12.cc index 5c75100d00..e82e3f9beb 100644 --- a/examples/step-12/step-12.cc +++ b/examples/step-12/step-12.cc @@ -115,6 +115,27 @@ namespace Step12 values[i]=0.; } } + + + // Finally, a function that computes and returns the wind field + // $\beta=\beta(\mathbf x)$. As explained in the introduction, we + // will use a rotational field around the origin in 2d. In 3d, we + // simply leave the $z$-component unset (i.e., at zero), whereas + // the function can not be used in 1d in its current implementation: + template + Tensor<1,dim> beta (const Point &p) + { + Assert (dim >= 2, ExcNotImplemented()); + + Point wind_field; + wind_field(0) = -p(1); + wind_field(1) = p(0); + wind_field /= wind_field.norm(); + + return wind_field; + } + + // @sect3{The AdvectionProblem class} // // After this preparations, we proceed with the main class of this program, @@ -367,22 +388,17 @@ namespace Step12 for (unsigned int point=0; point beta; - beta(0) = -fe_face_values.quadrature_point(point)(1); - beta(1) = fe_face_values.quadrature_point(point)(0); - beta /= beta.norm(); - - const double beta_n=beta * normals[point]; - if (beta_n>0) + const double beta_dot_n = beta(fe_face_values.quadrature_point(point)) * normals[point]; + if (beta_dot_n>0) for (unsigned int i=0; i beta; - beta(0) = -fe_face_values.quadrature_point(point)(1); - beta(1) = fe_face_values.quadrature_point(point)(0); - beta /= beta.norm(); - - const double beta_dot_n = beta * normals[point]; + const double beta_dot_n = beta(fe_face_values.quadrature_point(point)) * normals[point]; if (beta_dot_n>0) { // This term we've already seen: