diff options
Diffstat (limited to 'src/gr/expansion.cc')
| -rw-r--r-- | src/gr/expansion.cc | 491 |
1 files changed, 437 insertions, 54 deletions
diff --git a/src/gr/expansion.cc b/src/gr/expansion.cc index 5a03783..cdbd57c 100644 --- a/src/gr/expansion.cc +++ b/src/gr/expansion.cc @@ -37,9 +37,11 @@ #include <stdio.h> #include <assert.h> #include <math.h> +#include <string.h> #include "util_Table.h" #include "cctk.h" +#include "cctk_Arguments.h" #include "config.h" #include "stdc.h" @@ -47,6 +49,10 @@ #include "../jtutil/array.hh" #include "../jtutil/cpm_map.hh" #include "../jtutil/linear_map.hh" +using jtutil::error_exit; +using jtutil::pow2; +using jtutil::pow3; +using jtutil::pow4; #include "../patch/coords.hh" #include "../patch/grid.hh" @@ -65,9 +71,6 @@ // all the code in this file is inside this namespace namespace AHFinderDirect { -using jtutil::error_exit; -using jtutil::pow2; -using jtutil::pow4; //****************************************************************************** //****************************************************************************** @@ -99,8 +102,13 @@ bool geometry_is_finite(patch_system& ps, const struct error_info& error_info, bool initial_flag, bool print_msg_flag); -bool compute_Theta(patch_system& ps, fp add_to_expansion, - bool Jacobian_flag, jtutil::norm<fp>* Theta_norms_ptr, +bool compute_Theta(patch_system& ps, const struct what_to_compute& comput_info, + bool Jacobian_flag, + jtutil::norm<fp>* Theta_norms_ptr, + jtutil::norm<fp>* expansion_Theta_norms_ptr, + jtutil::norm<fp>* inner_expansion_Theta_norms_ptr, + jtutil::norm<fp>* product_expansion_Theta_norms_ptr, + jtutil::norm<fp>* mean_curvature_Theta_norms_ptr, const struct error_info& error_info, bool initial_flag, bool print_msg_flag); } @@ -172,13 +180,18 @@ bool compute_Theta(patch_system& ps, fp add_to_expansion, // succeeded or failed, and if the latter, what caused the failure. // enum expansion_status - expansion(patch_system* ps_ptr, fp add_to_expansion, + expansion(patch_system* ps_ptr, + const struct what_to_compute& compute_info, const struct cactus_grid_info& cgi, const struct geometry_info& gi, const struct error_info& error_info, bool initial_flag, bool Jacobian_flag /* = false */, bool print_msg_flag /* = false */, - jtutil::norm<fp>* Theta_norms_ptr /* = NULL */) + jtutil::norm<fp>* Theta_norms_ptr /* = NULL */, + jtutil::norm<fp>* expansion_Theta_norms_ptr /* = NULL */, + jtutil::norm<fp>* inner_expansion_Theta_norms_ptr /* = NULL */, + jtutil::norm<fp>* product_expansion_Theta_norms_ptr /* = NULL */, + jtutil::norm<fp>* mean_curvature_Theta_norms_ptr /* = NULL */) { const bool active_flag = (ps_ptr != NULL); if (print_msg_flag) @@ -230,6 +243,8 @@ if (gi.hardwire_Schwarzschild_EF_geometry) if (active_flag) then { + + if (gi.check_that_geometry_is_finite && !geometry_is_finite(*ps_ptr, error_info, initial_flag, @@ -237,15 +252,134 @@ if (active_flag) then return expansion_failure__geometry_nonfinite; // *** ERROR RETURN *** + // Ensure that there is a norm object + const bool want_norms = Theta_norms_ptr; + jtutil::norm<fp> norms; + if (compute_info.surface_selection != selection_definition) + then if (! Theta_norms_ptr) Theta_norms_ptr = &norms; + + // compute remaining gridfns --> $\Theta$ // and optionally also the Jacobian coefficients // by algebraic ops and angular finite differencing - if (!compute_Theta(*ps_ptr, add_to_expansion, + what_to_compute this_compute_info (compute_info); + this_compute_info.surface_selection = selection_definition; + if (!compute_Theta(*ps_ptr, this_compute_info, Jacobian_flag, Theta_norms_ptr, + expansion_Theta_norms_ptr, + inner_expansion_Theta_norms_ptr, + product_expansion_Theta_norms_ptr, + mean_curvature_Theta_norms_ptr, error_info, initial_flag, print_msg_flag)) then return expansion_failure__gij_not_positive_definite; // *** ERROR RETURN *** + + if (compute_info.surface_selection != selection_definition) { + // + // Apply correction to find a surface by its areal radius + // + // get mean expansion + fp mean_expansion; + fp areal_radius; + switch (compute_info.surface_selection) { + case selection_mean_coordinate_radius: { + const int np = ps_ptr->N_grid_points(); + fp sum_expansion = 0; + fp sum_radius = 0; + for (int pn = 0; pn < ps_ptr->N_patches(); ++pn) { + patch& p = ps_ptr->ith_patch(pn); + for (int irho = p.min_irho(); irho <= p.max_irho(); ++irho) { + for (int isigma = p.min_isigma(); isigma <= p.max_isigma(); ++isigma) { + sum_expansion += p.gridfn(gfns::gfn__Theta, irho,isigma); + sum_radius += p.ghosted_gridfn(gfns::gfn__h, irho,isigma); + } + } + } + mean_expansion = sum_expansion / np; + areal_radius = sum_radius / np; + break; + } + case selection_areal_radius: { + // get surface area + const fp area = ps_ptr->integrate_gridfn + (gfns::gfn__one, true, true, true, + gfns::gfn__h, + gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13, + gfns::gfn__g_dd_22, gfns::gfn__g_dd_23, + gfns::gfn__g_dd_33, + patch::integration_method__automatic_choice); + mean_expansion = Theta_norms_ptr->mean(); + areal_radius = sqrt(area / (4.0*PI)); + break; + } + case selection_expansion_mean_coordinate_radius: { + const int np = ps_ptr->N_grid_points(); + fp sum_expansion = 0; + fp sum_radius = 0; + for (int pn = 0; pn < ps_ptr->N_patches(); ++pn) { + patch& p = ps_ptr->ith_patch(pn); + for (int irho = p.min_irho(); irho <= p.max_irho(); ++irho) { + for (int isigma = p.min_isigma(); isigma <= p.max_isigma(); ++isigma) { + sum_expansion += p.gridfn(gfns::gfn__Theta, irho,isigma); + sum_radius += p.ghosted_gridfn(gfns::gfn__h, irho,isigma); + } + } + } + mean_expansion = sum_expansion / np; + areal_radius = mean_expansion * sum_radius / np; + break; + } + case selection_expansion_areal_radius: { + // get surface area + const fp area = ps_ptr->integrate_gridfn + (gfns::gfn__one, true, true, true, + gfns::gfn__h, + gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13, + gfns::gfn__g_dd_22, gfns::gfn__g_dd_23, + gfns::gfn__g_dd_33, + patch::integration_method__automatic_choice); + mean_expansion = Theta_norms_ptr->mean(); + areal_radius = mean_expansion * sqrt(area / (4.0*PI)); + break; + } + default: + assert (0); + } // switch areal_radius_definition + + if (! ps_ptr->N_additional_points()) { + // calculate correction + const fp correction + = (- mean_expansion + + areal_radius - compute_info.desired_value); + // apply correction + ps_ptr->add_to_gridfn(-correction, gfns::gfn__Theta); + } else { + const int np = ps_ptr->N_grid_points(); + const int gnp = ps_ptr->ghosted_N_grid_points(); + // apply correction + const fp correction + = ps_ptr->ghosted_gridfn_data(gfns::gfn__h)[gnp]; + ps_ptr->add_to_gridfn(-correction, gfns::gfn__Theta); + ps_ptr->gridfn_data(gfns::gfn__Theta)[np] + = (mean_expansion - correction + - areal_radius + compute_info.desired_value); + } + if (want_norms) { + // recalculate norms + ps_ptr->gridfn_norms (gfns::gfn__Theta, *Theta_norms_ptr); + } // if want_norms + } else { + // + // do not apply correction + // + if (ps_ptr->N_additional_points()) { + const int np = ps_ptr->N_grid_points(); + const int gnp = ps_ptr->ghosted_N_grid_points(); + ps_ptr->gridfn_data(gfns::gfn__Theta)[np] + = ps_ptr->ghosted_gridfn_data(gfns::gfn__h)[gnp]; + } + } // if use-areal-radius } return expansion_success; // *** NORMAL RETURN *** @@ -256,25 +390,19 @@ return expansion_success; // *** NORMAL RETURN *** //****************************************************************************** // -// This function sets up the xyz-position gridfns: -// * global_[xyz] (used by interplate_geometry() ) -// * global_{xx,xy,xz,yy,yz,zz} (used by higher-level code to compute -// quadrupole moments of horizons) -// -// Bugs: -// * We initialize the global_{xx,xy,xz,yy,yz,zz} gridfns every time -// this function is called, i.e. at each horizon-finding Newton -// iteration, even though they're only needed at the end of the -// horizon-finding process. In practice the extra cost is small, -// though, so it's probably not worth fixing this... +// This function sets up the global xyz positions of the grid points +// in the gridfns global_[xyz]. These will be used by interplate_geometry(). // namespace { -void setup_xyz_posns(patch_system& ps, bool print_msg_flag) +void setup_xyz_posns(patch_system& ps, const bool print_msg_flag) { if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, " xyz positions and derivative coefficients"); +#ifdef GEOMETRY_INTERP_DEBUG2 + printf ("AH exp\n"); +#endif for (int pn = 0 ; pn < ps.N_patches() ; ++pn) { patch& p = ps.ith_patch(pn); @@ -285,12 +413,15 @@ if (print_msg_flag) isigma <= p.max_isigma() ; ++isigma) { - const fp r = p.ghosted_gridfn(gfns::gfn__h, irho,isigma); + const fp r = p.ghosted_gridfn(gfns::gfn__h, irho,isigma); const fp rho = p.rho_of_irho(irho); const fp sigma = p.sigma_of_isigma(isigma); - fp local_x, local_y, local_z; p.xyz_of_r_rho_sigma(r,rho,sigma, local_x,local_y,local_z); +#ifdef GEOMETRY_INTERP_DEBUG2 + printf (" pn=%d irho=%d isigma=%d x=%g y=%g z=%g\n", + pn, irho, isigma, local_x, local_y, local_z); +#endif const fp global_x = ps.origin_x() + local_x; const fp global_y = ps.origin_y() + local_y; @@ -353,12 +484,15 @@ if (print_msg_flag) // global_[xyz] # xyz positions of grid points // // Inputs (Cactus 3-D gridfns): +// ahmask # excision mask // gxx,gxy,gxz,gyy,gyz,gzz # 3-metric $g_{ij}$ // # (may be either physical or conformal) // kxx,kxy,kxz,kyy,kyz,kzz # extrinsic curvature $K_{ij}$ // psi # optional conformal factor $\psi$ // // Outputs (angular gridfns, all on the nominal grid): +// mask # excision mask +// partial_d_mask_[123] # derivatives of the mask // g_dd_{11,12,13,22,23,33} # $\stored{g}_{ij}$ // K_dd_{11,12,13,22,23,33} # $K_{ij}$ // partial_d_g_dd_[123]{11,12,13,22,23,33} # $\partial_k \stored{g}_{ij}$ @@ -428,6 +562,7 @@ const void* const interp_coords[N_GRID_DIMS] const CCTK_INT input_array_variable_indices[] = { + cgi.mask_varindex, cgi.g_dd_11_varindex, cgi.g_dd_12_varindex, cgi.g_dd_13_varindex, cgi.g_dd_22_varindex, cgi.g_dd_23_varindex, cgi.g_dd_33_varindex, @@ -436,19 +571,6 @@ const CCTK_INT input_array_variable_indices[] cgi.K_dd_33_varindex, cgi.psi_varindex, }; -const CCTK_INT input_array_type_codes[] - = { - // g_ij - CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, - CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, - CCTK_VARIABLE_REAL, - // K_ij - CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, - CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, - CCTK_VARIABLE_REAL, - // psi - CCTK_VARIABLE_REAL, - }; const int N_input_arrays_for_psi = 1; const int N_input_arrays_dim = sizeof(input_array_variable_indices) / sizeof(input_array_variable_indices[0]); @@ -463,6 +585,8 @@ const int N_input_arrays_use const CCTK_INT output_array_type_codes[] = { + // mask $\partial_x$ mask $\partial_y$ mask $\partial_z$ mask + CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, // $g_{ij}$ $\partial_x g_{ij}$ $\partial_y g_{ij}$ $\partial_z g_{ij}$ CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, @@ -480,18 +604,20 @@ const CCTK_INT output_array_type_codes[] const CCTK_INT operand_indices[] = { - 0, 0, 0, 0, // g_dd_11, partial_[xyz] g_dd_11 - 1, 1, 1, 1, // g_dd_12, partial_[xyz] g_dd_12 - 2, 2, 2, 2, // g_dd_13, partial_[xyz] g_dd_13 - 3, 3, 3, 3, // g_dd_22, partial_[xyz] g_dd_22 - 4, 4, 4, 4, // g_dd_23, partial_[xyz] g_dd_23 - 5, 5, 5, 5, // g_dd_33, partial_[xyz] g_dd_33 - 6, 7, 8, 9, 10, 11, // K_dd_{11,12,13,22,23,33} - 12, 12, 12, 12, // psi, partial_[xyz] psi + 0, 0, 0, 0, // mask, partial_[xyz] mask + 1, 1, 1, 1, // g_dd_11, partial_[xyz] g_dd_11 + 2, 2, 2, 2, // g_dd_12, partial_[xyz] g_dd_12 + 3, 3, 3, 3, // g_dd_13, partial_[xyz] g_dd_13 + 4, 4, 4, 4, // g_dd_22, partial_[xyz] g_dd_22 + 5, 5, 5, 5, // g_dd_23, partial_[xyz] g_dd_23 + 6, 6, 6, 6, // g_dd_33, partial_[xyz] g_dd_33 + 7, 8, 9, 10, 11, 12, // K_dd_{11,12,13,22,23,33} + 13, 13, 13, 13, // psi, partial_[xyz] psi }; #define DERIV(x) x const CCTK_INT operation_codes[] = { + DERIV(0), DERIV(1), DERIV(2), DERIV(3), // mask, partial_[xyz] mask DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_11, partial_[xyz] g_dd_11 DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_12, partial_[xyz] g_dd_12 DERIV(0), DERIV(1), DERIV(2), DERIV(3), // g_dd_13, partial_[xyz] g_dd_13 @@ -505,6 +631,10 @@ const CCTK_INT operation_codes[] void* const output_arrays[] = { + CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__mask)), + CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_mask_1)), + CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_mask_2)), + CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_mask_3)), CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__g_dd_11)), CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_111)), CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_211)), @@ -625,6 +755,21 @@ if (print_msg_flag) " calling geometry interpolator (%s%d points)", (active_flag ? "" : "dummy: "), N_interp_points); +#ifdef GEOMETRY_INTERP_DEBUG2 + { + printf("AHFinderDirect:: proc %d: CCTK_InterpGridArrays() coordinates are:\n", + int(CCTK_MyProc(cgi.GH))); + for (int pt = 0 ; pt < N_interp_points ; ++ pt) + { + printf(" pt=%d [x,y,z]=[%g,%g,%g]\n", + pt, + double(((const CCTK_REAL*)interp_coords[0])[pt]), + double(((const CCTK_REAL*)interp_coords[1])[pt]), + double(((const CCTK_REAL*)interp_coords[2])[pt])); + } + } +#endif /* GEOMETRY_INTERP_DEBUG2 */ + #ifdef GEOMETRY_INTERP_DEBUG printf("AHFinderDirect:: proc %d: initializing interpolator outputs to 999.999\n", int(CCTK_MyProc(cgi.GH))); @@ -676,8 +821,11 @@ fflush(stdout); { for (int pt = 0 ; pt < N_interp_points ; pt = 2*pt + (pt == 0)) { - printf("AHFinderDirect:: proc %d: CCTK_InterpGridArrays() results for pt=%d:\n", - int(CCTK_MyProc(cgi.GH)), pt); + printf("AHFinderDirect:: proc %d: CCTK_InterpGridArrays() results for pt=%d at [x,y,z]=[%g,%g,%g]:\n", + int(CCTK_MyProc(cgi.GH)), pt, + double(((const CCTK_REAL*)interp_coords[0])[pt]), + double(((const CCTK_REAL*)interp_coords[1])[pt]), + double(((const CCTK_REAL*)interp_coords[2])[pt])); for (int out = 0 ; out < N_output_arrays_use ; ++out) { const CCTK_REAL* const out_ptr @@ -743,8 +891,8 @@ if (status == CCTK_ERROR_INTERP_POINT_OUTSIDE) CCTK_VWarn(warn_level, __LINE__, __FILE__, CCTK_THORNSTRING, "interpolate_geometry():\n" -" one or more points on the trial horizon surface point" -" is/are outside the grid (or too close to the grid boundary)" +" one or more points on the trial horizon surface point\n" +" is/are outside the grid (or too close to the grid boundary)\n" " (in a single-processor run, this may also mean that\n" " driver::ghost_size is too small for this geometry interpolator)\n"); } @@ -763,6 +911,51 @@ else if (status < 0) "***** interpolate_geometry(): error return %d from interpolator!\n", status); /*NOTREACHED*/ +// +// ***** check the interpolated excision mask ***** +// +{ +bool did_use_excised_gridpoint = false; +if (active_flag) + then { + for (int pn = 0 ; pn < ps_ptr->N_patches() ; ++pn) + { + patch& p = ps_ptr->ith_patch(pn); + + for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho) + for (int isigma = p.min_isigma() ; isigma <= p.max_isigma() ; ++isigma) + { + const fp m = p.gridfn(gfns::gfn__mask, irho,isigma); + const fp m1 = p.gridfn(gfns::gfn__partial_d_mask_1, irho,isigma); + const fp m2 = p.gridfn(gfns::gfn__partial_d_mask_2, irho,isigma); + const fp m3 = p.gridfn(gfns::gfn__partial_d_mask_3, irho,isigma); + if (fabs(m) > 1.0e-12 + || fabs(m1) > 1.0e-12 || fabs(m2) > 1.0e-12 || fabs(m3) > 1.0e-12) + then did_use_excised_gridpoint = true; + } + } + } +if (did_use_excised_gridpoint) + then { + if (print_msg_flag) + then { + // see if we can get further info + const int warn_level + = initial_flag + ? error_info.warn_level__point_outside__initial + : error_info.warn_level__point_outside__subsequent; + + + CCTK_VWarn(warn_level, __LINE__, __FILE__, CCTK_THORNSTRING, +"interpolate_geometry():\n" +" one or more points on the trial horizon surface point\n" +" is/are in an excised region (or too close to the excision boundary)\n"); + } + + return expansion_failure__surface_in_excised_region; // *** ERROR RETURN *** + } +} + return expansion_success; // *** NORMAL RETURN *** } } @@ -1166,14 +1359,59 @@ return true; // *** no check possible *** // g_ij isn't positive definite). // namespace { -bool compute_Theta(patch_system& ps, fp add_to_expansion, - bool Jacobian_flag, jtutil::norm<fp>* Theta_norms_ptr, +bool compute_Theta(patch_system& ps, const struct what_to_compute& compute_info, + bool Jacobian_flag, + jtutil::norm<fp>* Theta_norms_ptr, + jtutil::norm<fp>* expansion_Theta_norms_ptr, + jtutil::norm<fp>* inner_expansion_Theta_norms_ptr, + jtutil::norm<fp>* product_expansion_Theta_norms_ptr, + jtutil::norm<fp>* mean_curvature_Theta_norms_ptr, const struct error_info& error_info, bool initial_flag, bool print_msg_flag) { if (print_msg_flag) then CCTK_VInfo(CCTK_THORNSTRING, " computing Theta(h)"); +#if 0 + fp mean_radius, areal_radius; + switch (compute_info.surface_modification) { + case modification_none: + case modification_radius: + case modification_radius2: + // do nothing + break; + case modification_mean_radius: { + // get average coordinate radius + const int np = ps.N_grid_points(); + fp sum_radius = 0; + for (int pn = 0; pn < ps.N_patches(); ++pn) { + patch& p = ps.ith_patch(pn); + for (int irho = p.min_irho(); irho <= p.max_irho(); ++irho) { + for (int isigma = p.min_isigma(); isigma <= p.max_isigma(); ++isigma) { + sum_radius += p.ghosted_gridfn(gfns::gfn__h, irho,isigma); + } + } + } + mean_radius = sum_radius / np; + break; + } + case modification_areal_radius: { + // get surface area + const fp area = ps.integrate_gridfn + (gfns::gfn__one, true, true, true, + gfns::gfn__h, + gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13, + gfns::gfn__g_dd_22, gfns::gfn__g_dd_23, + gfns::gfn__g_dd_33, + patch::integration_method__automatic_choice); + areal_radius = sqrt(area / (4.0*PI)); + break; + } + default: + assert (0); + } +#endif + for (int pn = 0 ; pn < ps.N_patches() ; ++pn) { patch& p = ps.ith_patch(pn); @@ -1284,24 +1522,169 @@ CCTK_VWarn(warn_level, __LINE__, __FILE__, CCTK_THORNSTRING, return false; // *** ERROR RETURN *** } + assert (compute_info.surface_selection == selection_definition); + // compute H via equation (14) of my 1996 horizon finding paper const fp sqrt_Theta_D = sqrt(Theta_D); - Theta = + Theta_A/(Theta_D*sqrt_Theta_D) - + Theta_B/sqrt_Theta_D - + Theta_C/Theta_D - - K - + add_to_expansion; + const fp Theta_X = + Theta_A/(Theta_D*sqrt_Theta_D) + + Theta_B/sqrt_Theta_D; + const fp Theta_Y = + Theta_C/Theta_D + - K; + +#define mean_curvature p.gridfn(gfns::gfn__mean_curvature, irho,isigma) + mean_curvature = Theta_X; +#undef mean_curvature + + switch (compute_info.surface_definition) { + case definition_expansion: + Theta = + Theta_X + Theta_Y; + break; + case definition_inner_expansion: + Theta = - Theta_X + Theta_Y; + break; + case definition_mean_curvature: + Theta = + Theta_X; + break; + case definition_expansion_product: + Theta = (+ Theta_X + Theta_Y) * (- Theta_X + Theta_Y); + break; + default: + assert (0); + } + + switch (compute_info.surface_modification) { + case modification_none: + // do nothing + break; + case modification_radius: + // multiply by radius + Theta *= r; + break; + case modification_radius2: + // multiply by radius^2 + Theta *= pow2(r); + break; +#if 0 + case modification_mean_radius: + // multiply by average coordinate radius + Theta *= mean_radius; + break; + case modification_areal_radius: + // multiply by areal radius + Theta *= areal_radius; + break; +#endif + default: + assert (0); + } + + Theta -= compute_info.desired_value; + // update running norms of Theta(h) function if (Theta_norms_ptr != NULL) then Theta_norms_ptr->data(Theta); + if (expansion_Theta_norms_ptr != NULL) + then expansion_Theta_norms_ptr->data(+ Theta_X + Theta_Y); + + if (inner_expansion_Theta_norms_ptr != NULL) + then inner_expansion_Theta_norms_ptr->data(- Theta_X + Theta_Y); + + if (product_expansion_Theta_norms_ptr != NULL) + then product_expansion_Theta_norms_ptr->data((+ Theta_X + Theta_Y) * (- Theta_X + Theta_Y)); + + if (mean_curvature_Theta_norms_ptr != NULL) + then mean_curvature_Theta_norms_ptr->data(+ Theta_X); + + fp partial_Theta_X_wrt_partial_d_h_1; + fp partial_Theta_X_wrt_partial_d_h_2; + fp partial_Theta_X_wrt_partial_dd_h_11; + fp partial_Theta_X_wrt_partial_dd_h_12; + fp partial_Theta_X_wrt_partial_dd_h_22; + fp partial_Theta_Y_wrt_partial_d_h_1; + fp partial_Theta_Y_wrt_partial_d_h_2; + fp partial_Theta_Y_wrt_partial_dd_h_11; + fp partial_Theta_Y_wrt_partial_dd_h_12; + fp partial_Theta_Y_wrt_partial_dd_h_22; + if (Jacobian_flag) then { // partial_Theta_wrt_partial_d_h, // partial_Theta_wrt_partial_dd_h #include "../gr.cg/expansion_Jacobian.c" } + + if (Jacobian_flag) { + switch (compute_info.surface_definition) { + + case definition_expansion: + partial_Theta_wrt_partial_d_h_1 + = (+ partial_Theta_X_wrt_partial_d_h_1 + + partial_Theta_Y_wrt_partial_d_h_1); + partial_Theta_wrt_partial_d_h_2 + = (+ partial_Theta_X_wrt_partial_d_h_2 + + partial_Theta_Y_wrt_partial_d_h_2); + partial_Theta_wrt_partial_dd_h_11 + = (+ partial_Theta_X_wrt_partial_dd_h_11 + + partial_Theta_Y_wrt_partial_dd_h_11); + partial_Theta_wrt_partial_dd_h_12 + = (+ partial_Theta_X_wrt_partial_dd_h_12 + + partial_Theta_Y_wrt_partial_dd_h_12); + partial_Theta_wrt_partial_dd_h_22 + = (+ partial_Theta_X_wrt_partial_dd_h_22 + + partial_Theta_Y_wrt_partial_dd_h_22); + break; + + case definition_inner_expansion: + partial_Theta_wrt_partial_d_h_1 + = (- partial_Theta_X_wrt_partial_d_h_1 + + partial_Theta_Y_wrt_partial_d_h_1); + partial_Theta_wrt_partial_d_h_2 + = (- partial_Theta_X_wrt_partial_d_h_2 + + partial_Theta_Y_wrt_partial_d_h_2); + partial_Theta_wrt_partial_dd_h_11 + = (- partial_Theta_X_wrt_partial_dd_h_11 + + partial_Theta_Y_wrt_partial_dd_h_11); + partial_Theta_wrt_partial_dd_h_12 + = (- partial_Theta_X_wrt_partial_dd_h_12 + + partial_Theta_Y_wrt_partial_dd_h_12); + partial_Theta_wrt_partial_dd_h_22 + = (- partial_Theta_X_wrt_partial_dd_h_22 + + partial_Theta_Y_wrt_partial_dd_h_22); + break; + + case definition_mean_curvature: + partial_Theta_wrt_partial_d_h_1 + = + partial_Theta_X_wrt_partial_d_h_1; + partial_Theta_wrt_partial_d_h_2 + = + partial_Theta_X_wrt_partial_d_h_2; + partial_Theta_wrt_partial_dd_h_11 + = + partial_Theta_X_wrt_partial_dd_h_11; + partial_Theta_wrt_partial_dd_h_12 + = + partial_Theta_X_wrt_partial_dd_h_12; + partial_Theta_wrt_partial_dd_h_22 + = + partial_Theta_X_wrt_partial_dd_h_22; + break; + + case definition_expansion_product: { +#define f(x,y,dx,dy) (- x*x + y*y) +#define df(x,y,dx,dy) (- 2*x*dx + 2*y*dy) + partial_Theta_wrt_partial_d_h_1 = df(Theta_X, Theta_Y, partial_Theta_X_wrt_partial_d_h_1 , partial_Theta_Y_wrt_partial_d_h_1 ); + partial_Theta_wrt_partial_d_h_2 = df(Theta_X, Theta_Y, partial_Theta_X_wrt_partial_d_h_2 , partial_Theta_Y_wrt_partial_d_h_2 ); + partial_Theta_wrt_partial_dd_h_11 = df(Theta_X, Theta_Y, partial_Theta_X_wrt_partial_dd_h_11, partial_Theta_Y_wrt_partial_dd_h_11); + partial_Theta_wrt_partial_dd_h_12 = df(Theta_X, Theta_Y, partial_Theta_X_wrt_partial_dd_h_12, partial_Theta_Y_wrt_partial_dd_h_12); + partial_Theta_wrt_partial_dd_h_22 = df(Theta_X, Theta_Y, partial_Theta_X_wrt_partial_dd_h_22, partial_Theta_Y_wrt_partial_dd_h_22); +#undef f +#undef df + break; + } + + default: + assert (0); + } + } + } } } |