Skip to content

REFAC: extract stress/strain/rotation into C APIs and call them from Python #277

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Dengda98 merged 1 commit into from
Aug 4, 2026
Merged

REFAC: extract stress/strain/rotation into C APIs and call them from Python #277

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions pygrt/C_extension/include/grt.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,7 @@
#include "grt/dynamic/grn.h"
#include "grt/dynamic/grnspec.h"
#include "grt/dynamic/layer.h"
#include "grt/dynamic/dyn_postprocess.h"
#include "grt/dynamic/signals.h"
#include "grt/dynamic/source.h"
#include "grt/dynamic/syn.h"
Expand All @@ -59,6 +60,7 @@

#include "grt/static/static_grn.h"
#include "grt/static/static_layer.h"
#include "grt/static/static_postprocess.h"
#include "grt/static/static_source.h"


Expand Down
35 changes: 35 additions & 0 deletions pygrt/C_extension/include/grt/dynamic/dyn_postprocess.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,35 @@
/**
* @file dyn_postprocess.h
* @brief 动态位移偏导后处理（应力/应变/旋转）
* @date 2026-07
*/

#pragma once

#include "grt/common/const.h"

/**
* 由动态位移偏导合成应变张量。
* 数组布局：u[分量][采样点]、upar[偏导方向][分量][采样点]、
* res[第二分量][第一分量][采样点]。
*/
void grt_compute_strain(
size_t npts, float dist, float *const u[GRT_CHANNEL_NUM],
float *const upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM],
float *const res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], bool rot2ZNE);

/** 由动态位移偏导合成旋转张量。数组布局同 grt_compute_strain()。 */
void grt_compute_rotation(
size_t npts, float dist, float *const u[GRT_CHANNEL_NUM],
float *const upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM],
float *const res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], bool rot2ZNE);

/**
* 在频域由动态位移偏导合成应力张量。
* 介质参数约定与命令行 stress 子模块的 SAC 头段一致。
*/
void grt_compute_stress(
size_t npts, float dt, float dist, float va, float vb, float rho,
float Qainv, float Qbinv, float *const u[GRT_CHANNEL_NUM],
float *const upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM],
float *const res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], bool rot2ZNE);
42 changes: 42 additions & 0 deletions pygrt/C_extension/include/grt/static/static_postprocess.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,42 @@
/**
* @file static_postprocess.h
* @brief 静态位移偏导后处理（应力/应变/旋转）
* @date 2026-07
*/

#pragma once

#include "grt/common/const.h"

/**
* 由静态位移偏导合成对称应力张量。
*
* 数组布局：u[分量][点]、upar[偏导方向][分量][点]、
* res[第二分量][第一分量][点]。仅写入 res 的上三角分量。
*/
void grt_static_compute_stress(
size_t nx, size_t ny, const real_t *xs, const real_t *ys,
real_t *const u[GRT_CHANNEL_NUM],
real_t *const upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM],
real_t *const res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM],
bool rot2ZNE, real_t mu, real_t lam);

/**
* 由静态位移偏导合成对称应变张量。
* 数组布局同 grt_static_compute_stress()。
*/
void grt_static_compute_strain(
size_t nx, size_t ny, const real_t *xs, const real_t *ys,
real_t *const u[GRT_CHANNEL_NUM],
real_t *const upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM],
real_t *const res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], bool rot2ZNE);

/**
* 由静态位移偏导合成反对称旋转张量。
* 数组布局同 grt_static_compute_stress()。仅写入 res 的非对角分量。
*/
void grt_static_compute_rotation(
size_t nx, size_t ny, const real_t *xs, const real_t *ys,
real_t *const u[GRT_CHANNEL_NUM],
real_t *const upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM],
real_t *const res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], bool rot2ZNE);
82 changes: 52 additions & 30 deletions pygrt/C_extension/src/dynamic/grt_rotation.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,29 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){
GRTCheckOptionSet(argc > 1);
}

void grt_compute_rotation(
size_t npts, float dist, float *const u[GRT_CHANNEL_NUM],
float *const upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM],
float *const res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], bool rot2ZNE)
{
const char *chs = rot2ZNE ? GRT_ZNE_CODES : GRT_ZRT_CODES;

for(size_t i=0; i<npts; ++i){
// ZRT 联络项 u_θ/r，1e-5: km→cm
float ut_over_r = u[2][i] / dist * 1e-5f;

for(int c=0; c<GRT_CHANNEL_NUM; ++c){
for(int c2=c+1; c2<GRT_CHANNEL_NUM; ++c2){
float val = 0.5f * (upar[c2][c][i] - upar[c][c2][i]);
if(chs[c]=='R' && chs[c2]=='T'){
val -= 0.5f * ut_over_r;
}
res[c2][c][i] = val;
}
}
}
}


/** 子模块主函数 */
int rotation_main(int argc, char **argv){
Expand Down Expand Up @@ -90,46 +113,45 @@ int rotation_main(int argc, char **argv){
SACTRACE *outsac = grt_copy_SACTRACE(insac, true);
grt_free_SACTRACE(insac);

// ----------------------------------------------------------------------------------
// 循环3个分量
float *u[GRT_CHANNEL_NUM];
float *upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM];
float *res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM];
for(int c=0; c<GRT_CHANNEL_NUM; ++c){
GRT_SAFE_ASPRINTF(&s_filepath, "%s/%c.sac", Ctrl->s_synpath, chs[c]);
insac = grt_read_SACTRACE(s_filepath, false);
u[c] = insac->data;
insac->data = NULL;
grt_free_SACTRACE(insac);
for(int c2=0; c2<GRT_CHANNEL_NUM; ++c2){
GRT_SAFE_ASPRINTF(&s_filepath, "%s/%c%c.sac", Ctrl->s_synpath, tolower(chs[c2]), chs[c]);
insac = grt_read_SACTRACE(s_filepath, false);
upar[c2][c] = insac->data;
insac->data = NULL;
grt_free_SACTRACE(insac);
res[c2][c] = calloc(npts, sizeof(*res[c2][c]));
}
}
grt_compute_rotation(npts, dist, u, upar, res, rot2ZNE);

// 写出3个分量
for(int i1=0; i1<2; ++i1){
c1 = chs[i1];
for(int i2=i1+1; i2<3; ++i2){
c2 = chs[i2];

// 读取数据 u_{i,j}
GRT_SAFE_ASPRINTF(&s_filepath, "%s/%c%c.sac", Ctrl->s_synpath, tolower(c2), c1);
insac = grt_read_SACTRACE(s_filepath, false);

// 累加
for(int i=0; i<npts; ++i) outsac->data[i] += insac->data[i];

// 读取数据 u_{j,i}
GRT_SAFE_ASPRINTF(&s_filepath, "%s/%c%c.sac", Ctrl->s_synpath, tolower(c1), c2);
insac = grt_read_SACTRACE(s_filepath, false);

// 累加
for(int i=0; i<npts; ++i) outsac->data[i] = (outsac->data[i] - insac->data[i]) * 0.5f;

// 特殊情况需加上协变导数，1e-5是因为km->cm
if(c1=='R' && c2=='T'){
// 读取数据 u_T
GRT_SAFE_ASPRINTF(&s_filepath, "%s/T.sac", Ctrl->s_synpath);
insac = grt_read_SACTRACE(s_filepath, false);
for(int i=0; i<npts; ++i) outsac->data[i] -= 0.5f * insac->data[i] / dist * 1e-5;
}

// 保存到SAC
memcpy(outsac->data, res[i2][i1], sizeof(*outsac->data)*npts);
sprintf(outsac->hd.kcmpnm, "%c%c", c1, c2);
GRT_SAFE_ASPRINTF(&s_filepath, "%s/rotation_%c%c.sac", Ctrl->s_synpath, c1, c2);
grt_write_SACTRACE(s_filepath, outsac);

// 置零
for(int i=0; i<npts; ++i) outsac->data[i] = 0.0f;
}
}

grt_free_SACTRACE(insac);
for(int c=0; c<GRT_CHANNEL_NUM; ++c){
free(u[c]);
for(int c2=0; c2<GRT_CHANNEL_NUM; ++c2){
free(upar[c2][c]);
free(res[c2][c]);
}
}
grt_free_SACTRACE(outsac);
GRT_SAFE_FREE_PTR(s_filepath);

Expand Down
92 changes: 56 additions & 36 deletions pygrt/C_extension/src/dynamic/grt_strain.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,33 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){
GRTCheckOptionSet(argc > 1);
}

void grt_compute_strain(
size_t npts, float dist, float *const u[GRT_CHANNEL_NUM],
float *const upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM],
float *const res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], bool rot2ZNE)
{
const char *chs = rot2ZNE ? GRT_ZNE_CODES : GRT_ZRT_CODES;

for(size_t i=0; i<npts; ++i){
// ZRT 联络项 u/r，1e-5: km→cm
float ur_over_r = u[1][i] / dist * 1e-5f;
float ut_over_r = u[2][i] / dist * 1e-5f;

for(int c=0; c<GRT_CHANNEL_NUM; ++c){
for(int c2=c; c2<GRT_CHANNEL_NUM; ++c2){
float val = 0.5f * (upar[c2][c][i] + upar[c][c2][i]);
if(chs[c]=='R' && chs[c2]=='T'){
val -= 0.5f * ut_over_r;
}
else if(chs[c]=='T' && chs[c2]=='T'){
val += ur_over_r;
}
res[c2][c][i] = val;
}
}
}
}



/** 子模块主函数 */
Expand Down Expand Up @@ -88,52 +115,45 @@ int strain_main(int argc, char **argv){
SACTRACE *outsac = grt_copy_SACTRACE(insac, true);
grt_free_SACTRACE(insac);

// ----------------------------------------------------------------------------------
// 循环6个分量
float *u[GRT_CHANNEL_NUM];
float *upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM];
float *res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM];
for(int c=0; c<GRT_CHANNEL_NUM; ++c){
GRT_SAFE_ASPRINTF(&s_filepath, "%s/%c.sac", Ctrl->s_synpath, chs[c]);
insac = grt_read_SACTRACE(s_filepath, false);
u[c] = insac->data;
insac->data = NULL;
grt_free_SACTRACE(insac);
for(int c2=0; c2<GRT_CHANNEL_NUM; ++c2){
GRT_SAFE_ASPRINTF(&s_filepath, "%s/%c%c.sac", Ctrl->s_synpath, tolower(chs[c2]), chs[c]);
insac = grt_read_SACTRACE(s_filepath, false);
upar[c2][c] = insac->data;
insac->data = NULL;
grt_free_SACTRACE(insac);
res[c2][c] = calloc(npts, sizeof(*res[c2][c]));
}
}
grt_compute_strain(npts, dist, u, upar, res, rot2ZNE);

// 写出6个分量
for(int i1=0; i1<3; ++i1){
c1 = chs[i1];
for(int i2=i1; i2<3; ++i2){
c2 = chs[i2];

// 读取数据 u_{i,j}
GRT_SAFE_ASPRINTF(&s_filepath, "%s/%c%c.sac", Ctrl->s_synpath, tolower(c2), c1);
insac = grt_read_SACTRACE(s_filepath, false);

// 累加
for(int i=0; i<npts; ++i) outsac->data[i] += insac->data[i];

// 读取数据 u_{j,i}
GRT_SAFE_ASPRINTF(&s_filepath, "%s/%c%c.sac", Ctrl->s_synpath, tolower(c1), c2);
insac = grt_read_SACTRACE(s_filepath, false);

// 累加
for(int i=0; i<npts; ++i) outsac->data[i] = (outsac->data[i] + insac->data[i]) * 0.5f;

// 特殊情况需加上协变导数，1e-5是因为km->cm
if(c1=='R' && c2=='T'){
// 读取数据 u_T
GRT_SAFE_ASPRINTF(&s_filepath, "%s/T.sac", Ctrl->s_synpath);
insac = grt_read_SACTRACE(s_filepath, false);
for(int i=0; i<npts; ++i) outsac->data[i] -= 0.5f * insac->data[i] / dist * 1e-5;
}
else if(c1=='T' && c2=='T'){
// 读取数据 u_R
GRT_SAFE_ASPRINTF(&s_filepath, "%s/R.sac", Ctrl->s_synpath);
insac = grt_read_SACTRACE(s_filepath, false);
for(int i=0; i<npts; ++i) outsac->data[i] += insac->data[i] / dist * 1e-5;
}

// 保存到SAC
memcpy(outsac->data, res[i2][i1], sizeof(*outsac->data)*npts);
sprintf(outsac->hd.kcmpnm, "%c%c", c1, c2);
GRT_SAFE_ASPRINTF(&s_filepath, "%s/strain_%c%c.sac", Ctrl->s_synpath, c1, c2);
grt_write_SACTRACE(s_filepath, outsac);

// 置零
for(int i=0; i<npts; ++i) outsac->data[i] = 0.0f;
}
}

grt_free_SACTRACE(insac);
for(int c=0; c<GRT_CHANNEL_NUM; ++c){
free(u[c]);
for(int c2=0; c2<GRT_CHANNEL_NUM; ++c2){
free(upar[c2][c]);
free(res[c2][c]);
}
}
grt_free_SACTRACE(outsac);
GRT_SAFE_FREE_PTR(s_filepath);

Expand Down
Loading
Loading