CGH0S7
4a413d669d
New files provide C equivalents of Fortran diff_new_sh.f90 and kodiss_sh.f90: - fderivs_sh_c.C: first derivatives in shell (rho, sigma, R) coords - fdderivs_sh_c.C: second derivatives in shell coords - fderivs_shc_c.C: shell first derivs + chain rule to Cartesian - fdderivs_shc_c.C: shell second derivs + chain rule to Cartesian - kodiss_sh_c.C: Kreiss-Oliger dissipation on shell patches Also add symmetry_stbd() C implementation and shell fh indexing to share_func.h. Controlled by USE_CXX_SHELL_KERNELS switch (default 0, experimental). Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
116 lines
4.9 KiB
C
116 lines
4.9 KiB
C
#include "macrodef.h"
|
|
#include "share_func.h"
|
|
#include <cstddef>
|
|
|
|
/* Forward declarations */
|
|
void fderivs_sh(const int ex[3], const double *f,
|
|
double *fx, double *fy, double *fz,
|
|
const double *X, const double *Y, const double *Z,
|
|
double SYM1, double SYM2, double SYM3,
|
|
int Symmetry, int onoff, int sst);
|
|
|
|
void fdderivs_sh(const int ex[3], const double *f,
|
|
double *fxx, double *fxy, double *fxz,
|
|
double *fyy, double *fyz, double *fzz,
|
|
const double *X, const double *Y, const double *Z,
|
|
double SYM1, double SYM2, double SYM3,
|
|
int Symmetry, int onoff, int sst);
|
|
|
|
/*
|
|
* fdderivs_shc — shell second derivatives converted to Cartesian via chain rule.
|
|
*
|
|
* Calls fderivs_sh and fdderivs_sh internally, then applies:
|
|
* f_{,ij} = sum_{a,b} (dx^a/dx^i)(dx^b/dx^j) f_{,ab} + sum_a (d^2 x^a/dx^i dx^j) f_{,a}
|
|
* where a,b ∈ {rho, sigma, R}.
|
|
*/
|
|
extern "C" {
|
|
|
|
void f_fdderivs_shc(int *ex,
|
|
double *f,
|
|
double *fxx, double *fxy, double *fxz,
|
|
double *fyy, double *fyz, double *fzz,
|
|
double *crho, double *sigma, double *R,
|
|
double &SYM1, double &SYM2, double &SYM3,
|
|
int &Symmetry, int &Lev, int &sst,
|
|
double *drhodx, double *drhody, double *drhodz,
|
|
double *dsigmadx, double *dsigmady, double *dsigmadz,
|
|
double *dRdx, double *dRdy, double *dRdz,
|
|
double *drhodxx, double *drhodxy, double *drhodxz,
|
|
double *drhodyy, double *drhodyz, double *drhodzz,
|
|
double *dsigmadxx, double *dsigmadxy, double *dsigmadxz,
|
|
double *dsigmadyy, double *dsigmadyz, double *dsigmadzz,
|
|
double *dRdxx, double *dRdxy, double *dRdxz,
|
|
double *dRdyy, double *dRdyz, double *dRdzz)
|
|
{
|
|
const int ex3[3] = { ex[0], ex[1], ex[2] };
|
|
const size_t n = (size_t)ex[0] * (size_t)ex[1] * (size_t)ex[2];
|
|
|
|
double *gx = (double*)malloc(n * sizeof(double));
|
|
double *gy = (double*)malloc(n * sizeof(double));
|
|
double *gz = (double*)malloc(n * sizeof(double));
|
|
double *gxx = (double*)malloc(n * sizeof(double));
|
|
double *gxy = (double*)malloc(n * sizeof(double));
|
|
double *gxz = (double*)malloc(n * sizeof(double));
|
|
double *gyy = (double*)malloc(n * sizeof(double));
|
|
double *gyz = (double*)malloc(n * sizeof(double));
|
|
double *gzz = (double*)malloc(n * sizeof(double));
|
|
|
|
if (!gx||!gy||!gz||!gxx||!gxy||!gxz||!gyy||!gyz||!gzz) {
|
|
free(gx);free(gy);free(gz);free(gxx);free(gxy);free(gxz);free(gyy);free(gyz);free(gzz);
|
|
return;
|
|
}
|
|
|
|
fderivs_sh(ex3, f, gx, gy, gz, crho, sigma, R, SYM1, SYM2, SYM3, Symmetry, Lev, sst);
|
|
fdderivs_sh(ex3, f, gxx, gxy, gxz, gyy, gyz, gzz, crho, sigma, R, SYM1, SYM2, SYM3, Symmetry, Lev, sst);
|
|
|
|
for (size_t i = 0; i < n; ++i) {
|
|
const double rx=drhodx[i], ry=drhody[i], rz=drhodz[i];
|
|
const double sx=dsigmadx[i], sy=dsigmady[i], sz=dsigmadz[i];
|
|
const double Rx=dRdx[i], Ry=dRdy[i], Rz=dRdz[i];
|
|
const double rxx=drhodxx[i], rxy=drhodxy[i], rxz=drhodxz[i];
|
|
const double ryy=drhodyy[i], ryz=drhodyz[i], rzz=drhodzz[i];
|
|
const double sxx=dsigmadxx[i], sxy=dsigmadxy[i], sxz=dsigmadxz[i];
|
|
const double syy=dsigmadyy[i], syz=dsigmadyz[i], szz=dsigmadzz[i];
|
|
const double Rxx=dRdxx[i], Rxy=dRdxy[i], Rxz=dRdxz[i];
|
|
const double Ryy=dRdyy[i], Ryz=dRdyz[i], Rzz=dRdzz[i];
|
|
|
|
const double Gr=gx[i], Gs=gy[i], GR=gz[i];
|
|
const double Grr=gxx[i], Grs=gxy[i], GrR=gxz[i];
|
|
const double Gss=gyy[i], GsR=gyz[i], GRR=gzz[i];
|
|
|
|
/* fxx */
|
|
fxx[i] = rx*rx*Grr + sx*sx*Gss + Rx*Rx*GRR
|
|
+ 2.0*(rx*sx*Grs + rx*Rx*GrR + sx*Rx*GsR)
|
|
+ rxx*Gr + sxx*Gs + Rxx*GR;
|
|
|
|
/* fxy */
|
|
fxy[i] = rx*ry*Grr + sx*sy*Gss + Rx*Ry*GRR
|
|
+ rx*sy*Grs + ry*sx*Grs + rx*Ry*GrR + ry*Rx*GrR + sx*Ry*GsR + sy*Rx*GsR
|
|
+ rxy*Gr + sxy*Gs + Rxy*GR;
|
|
|
|
/* fxz */
|
|
fxz[i] = rx*rz*Grr + sx*sz*Gss + Rx*Rz*GRR
|
|
+ rx*sz*Grs + rz*sx*Grs + rx*Rz*GrR + rz*Rx*GrR + sx*Rz*GsR + sz*Rx*GsR
|
|
+ rxz*Gr + sxz*Gs + Rxz*GR;
|
|
|
|
/* fyy */
|
|
fyy[i] = ry*ry*Grr + sy*sy*Gss + Ry*Ry*GRR
|
|
+ 2.0*(ry*sy*Grs + ry*Ry*GrR + sy*Ry*GsR)
|
|
+ ryy*Gr + syy*Gs + Ryy*GR;
|
|
|
|
/* fyz */
|
|
fyz[i] = ry*rz*Grr + sy*sz*Gss + Ry*Rz*GRR
|
|
+ ry*sz*Grs + rz*sy*Grs + ry*Rz*GrR + rz*Ry*GrR + sy*Rz*GsR + sz*Ry*GsR
|
|
+ ryz*Gr + syz*Gs + Ryz*GR;
|
|
|
|
/* fzz */
|
|
fzz[i] = rz*rz*Grr + sz*sz*Gss + Rz*Rz*GRR
|
|
+ 2.0*(rz*sz*Grs + rz*Rz*GrR + sz*Rz*GsR)
|
|
+ rzz*Gr + szz*Gs + Rzz*GR;
|
|
}
|
|
|
|
free(gx);free(gy);free(gz);free(gxx);free(gxy);free(gxz);free(gyy);free(gyz);free(gzz);
|
|
}
|
|
|
|
} // extern "C"
|