197 constexpr int dim = 3;
198 const Point<dim> p = p_ - coordinate_system_offset;
199 const std::array<double, dim> sp =
201 const std::array<double, dim> sg = sgradient(sp, component);
204 const double cos_theta =
std::cos(sp[1]);
205 const double sin_theta =
std::sin(sp[1]);
206 const double cos_phi =
std::cos(sp[2]);
207 const double sin_phi =
std::sin(sp[2]);
211 cos_theta, sin_theta, cos_phi, sin_phi, unit_r, unit_theta, unit_phi);
217 res += unit_r * sg[0];
220 if (sg[1] * sin_phi != 0.)
223 res += unit_theta * sg[1] / (sp[0] * sin_phi);
229 res += unit_phi * sg[2] / sp[0];
253 constexpr int dim = 3;
254 const Point<dim> p = p_ - coordinate_system_offset;
255 const std::array<double, dim> sp =
257 const std::array<double, dim> sg = sgradient(sp, component);
258 const std::array<double, 6> sh = shessian(sp, component);
261 const double cos_theta =
std::cos(sp[1]);
262 const double sin_theta =
std::sin(sp[1]);
263 const double cos_phi =
std::cos(sp[2]);
264 const double sin_phi =
std::sin(sp[2]);
265 const double r = sp[0];
269 cos_theta, sin_theta, cos_phi, sin_phi, unit_r, unit_theta, unit_phi);
271 const double sin_phi2 = sin_phi * sin_phi;
272 const double r2 = r * r;
275 const double c_utheta2 =
276 sg[0] / r + ((sin_phi != 0.) ? (cos_phi * sg[2]) / (r2 * sin_phi) +
277 sh[1] / (r2 * sin_phi2) :
279 const double c_utheta_ur =
280 ((sin_phi != 0.) ? (r * sh[3] - sg[1]) / (r2 * sin_phi) : 0.);
281 const double c_utheta_uphi =
282 ((sin_phi != 0.) ? (sh[5] * sin_phi - cos_phi * sg[1]) / (r2 * sin_phi2) :
284 const double c_ur2 = sh[0];
285 const double c_ur_uphi = (r * sh[4] - sg[2]) / r2;
286 const double c_uphi2 = (sh[2] + r * sg[0]) / r2;
291 add_outer_product(res, c_utheta2, unit_theta);
293 add_outer_product(res, c_utheta_ur, unit_theta, unit_r);
295 add_outer_product(res, c_utheta_uphi, unit_theta, unit_phi);
297 add_outer_product(res, c_ur2, unit_r);
299 add_outer_product(res, c_ur_uphi, unit_r, unit_phi);
301 add_outer_product(res, c_uphi2, unit_phi);