From 3a950881538f908b4bcf94c69260ef1cb71e12ed Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 2 Mar 2018 10:12:44 -0700 Subject: [PATCH] Merge the two places where we compute beta in step-12 into one function. --- examples/step-12/step-12.cc | 41 +++++++++++++++++++++++-------------- 1 file changed, 26 insertions(+), 15 deletions(-) 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: -- 2.39.5