dev_yi #1
BIN
GCTL_logo.jpg
Normal file
BIN
GCTL_logo.jpg
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 33 KiB |
105
README.md
105
README.md
@ -1,2 +1,105 @@
|
||||
# gctl_potential
|
||||
|
||||
|
||||
# GCTL Potential
|
||||
|
||||
GCTL Potential 是地球物理计算工具与库（Geophysical Computational Tools & Library, GCTL）的重力和磁力场正演模拟模块。该模块提供了一系列用于计算重力场和磁场的核心功能。
|
||||
|
||||
## 功能特性
|
||||
|
||||
### 重力场计算 (G-Kernel)
|
||||
- 2D模型计算
|
||||
- `gkernel_triangle2d`: 二维三角形体
|
||||
- `gkernel_rectangle2d`: 二维矩形体
|
||||
- `gkernel_polygon`: 二维多边形体
|
||||
|
||||
- 3D模型计算
|
||||
- `gkernel_triangle`: 三维三角形体
|
||||
- `gkernel_sphere`: 球体
|
||||
- `gkernel_block`: 长方体
|
||||
- `gkernel_tetrahedron`: 四面体
|
||||
- `gkernel_tricone`: 三锥体
|
||||
- `gkernel_tesseroid`: 曲面格网体（需要 GCTL_POTENTIAL_TESS 选项）
|
||||
|
||||
### 磁场计算 (M-Kernel)
|
||||
- `mkernel_dipole`: 偶极子
|
||||
- `mkernel_block`: 长方体
|
||||
- `mkernel_triangle`: 三角形体
|
||||
- `mkernel_tricone`: 三锥体
|
||||
- `mkernel_tetrahedron`: 四面体
|
||||
- 标准实现
|
||||
- Ren2017 实现方法
|
||||
- `mkernel_tesseroid`: 曲面格网体（需要 GCTL_POTENTIAL_MAGTESS 选项）
|
||||
|
||||
### 网格与数据结构
|
||||
- `gm_data`: 基础数据结构
|
||||
- `gm_regular_grid`: 规则网格
|
||||
- `gm_regular_mesh_2d`: 二维规则网格
|
||||
- `gm_regular_mesh_3d`: 三维规则网格
|
||||
- `gm_regular_mesh_sph_3d`: 三维球面规则网格
|
||||
- `gm_tet_mesh`: 四面体网格
|
||||
|
||||
## 依赖
|
||||
- Eigen (矩阵运算库)
|
||||
- autodiff (自动微分库，用于梯度计算)
|
||||
|
||||
## 许可证
|
||||
GCTL Potential 采用双重许可证方案：
|
||||
1. GNU Lesser General Public License v2 或更高版本
|
||||
2. 商业许可证（需单独购买）
|
||||
|
||||
## 作者
|
||||
- Yi Zhang (yizhang-geo@zju.edu.cn)
|
||||
|
||||
## 安装
|
||||
使用 CMake 构建系统：
|
||||
```bash
|
||||
mkdir build
|
||||
cd build
|
||||
cmake ..
|
||||
make
|
||||
make install
|
||||
```
|
||||
|
||||
## 示例
|
||||
`example` 目录下提供了丰富的示例代码，展示了各种计算功能的使用方法：
|
||||
|
||||
### 重力场计算示例
|
||||
- `gobser_tri2d_ex.cpp`: 二维三角形体重力场计算
|
||||
- 演示了如何读取 GMSH 格式的网格文件
|
||||
- 计算重力位、重力分量和重力梯度张量
|
||||
- 结果输出为 DSV 格式
|
||||
- `gobser_tri2d_sph_ex.cpp`: 球面坐标系下的二维三角形体重力场计算
|
||||
|
||||
### 磁场计算示例
|
||||
- `mobser_tetra_ex.cpp`: 四面体磁场计算
|
||||
- 展示了标准实现和 Ren2017 方法的对比
|
||||
- 包含磁位、磁场分量和 ΔT 异常计算
|
||||
- 支持 NetCDF 格式输出
|
||||
- `mobser_tetra_ex2.cpp`: 四面体磁场计算的另一个示例
|
||||
- `mobser_block_ex.cpp`: 长方体磁场计算
|
||||
- `mobser_tricone_ex.cpp`: 三锥体磁场计算
|
||||
- `mobser_tesseroid_ex.cpp`: 曲面格网体磁场计算
|
||||
- `mobser_dipole_ex.cpp`: 偶极子磁场计算
|
||||
- `mobser_tri_ex.cpp`: 三角形体磁场计算
|
||||
- `mobser_tri_sph_ex.cpp`: 球面坐标系下的三角形体磁场计算
|
||||
|
||||
### 特殊功能示例
|
||||
- `power_spectrum_ex.cpp`: 功率谱计算示例
|
||||
- `read_IGRF_ex.cpp`: IGRF 地磁场模型数据读取示例
|
||||
- `mobser_block_gradient_ex.cpp`: 长方体磁场梯度计算
|
||||
|
||||
### 数据格式支持
|
||||
- 网格文件格式：
|
||||
- GMSH 格式 (.msh)
|
||||
- Tetgen 格式 (.node, .ele)
|
||||
- 输出格式：
|
||||
- DSV (Delimiter-Separated Values)
|
||||
- NetCDF
|
||||
- 文本格式
|
||||
|
||||
所有示例代码都包含详细的注释，并提供了完整的数据输入输出流程。用户可以根据需要修改示例代码中的参数来适应自己的应用场景。
|
||||
|
||||
## 注意事项
|
||||
- 某些功能（如 tesseroid 计算）需要在编译时启用特定选项
|
||||
- 商业用途请联系作者获取商业许可证
|
||||
|
||||
|
||||
10
data/cylinder2d/cylinder2d.geo
Normal file
10
data/cylinder2d/cylinder2d.geo
Normal file
@ -0,0 +1,10 @@
|
||||
|
||||
//+
|
||||
SetFactory("OpenCASCADE");
|
||||
Circle(1) = {0, -100, 0, 20, 0, 2*Pi};
|
||||
//+
|
||||
Curve Loop(1) = {1};
|
||||
//+
|
||||
Plane Surface(1) = {1};
|
||||
//+
|
||||
Physical Surface("Body", 2) = {1};
|
||||
202
data/cylinder2d/gravity_2d.txt
Normal file
202
data/cylinder2d/gravity_2d.txt
Normal file
@ -0,0 +1,202 @@
|
||||
X Y Vx Vz Vxx Vxz Vzx Vzz
|
||||
-500 0 0.0318575695325 0.00637151390649 5.88139745215e-05 2.45058227173e-05 2.45058227173e-05 -5.88139745215e-05
|
||||
-495 0 0.0321542531032 0.00649580870773 5.98638345862e-05 2.5216549831e-05 2.5216549831e-05 -5.98638345862e-05
|
||||
-490 0 0.0324562523665 0.00662372497276 6.09403884939e-05 2.59546200452e-05 2.59546200452e-05 -6.09403884939e-05
|
||||
-485 0 0.0327637028692 0.00675540265344 6.20444719184e-05 2.67213398872e-05 2.67213398872e-05 -6.20444719184e-05
|
||||
-480 0 0.0330767444064 0.006890988418 6.31769487241e-05 2.75180901884e-05 2.75180901884e-05 -6.31769487241e-05
|
||||
-475 0 0.0333955211651 0.00703063603475 6.43387117239e-05 2.83463309624e-05 2.83463309624e-05 -6.43387117239e-05
|
||||
-470 0 0.0337201818698 0.00717450678081 6.55306834159e-05 2.92076066434e-05 2.92076066434e-05 -6.55306834159e-05
|
||||
-465 0 0.0340508799335 0.00732276987816 6.67538166924e-05 3.01035517149e-05 3.01035517149e-05 -6.67538166924e-05
|
||||
-460 0 0.0343877736109 0.00747560295888 6.80090955104e-05 3.10358967607e-05 3.10358967607e-05 -6.80090955104e-05
|
||||
-455 0 0.0347310261554 0.00763319256163 6.92975355123e-05 3.20064749733e-05 3.20064749733e-05 -6.92975355123e-05
|
||||
-450 0 0.0350808059793 0.00779573466206 7.06201845857e-05 3.3017229157e-05 3.3017229157e-05 -7.06201845857e-05
|
||||
-445 0 0.0354372868156 0.00796343523946 7.19781233457e-05 3.40702192675e-05 3.40702192675e-05 -7.19781233457e-05
|
||||
-440 0 0.0358006478832 0.00813651088255 7.33724655224e-05 3.51676305336e-05 3.51676305336e-05 -7.33724655224e-05
|
||||
-435 0 0.0361710740532 0.00831518943751 7.48043582351e-05 3.63117822092e-05 3.63117822092e-05 -7.48043582351e-05
|
||||
-430 0 0.0365487560157 0.00849971070132 7.62749821273e-05 3.75051370094e-05 3.75051370094e-05 -7.62749821273e-05
|
||||
-425 0 0.0369338904481 0.00869032716427 7.77855513392e-05 3.87503112898e-05 3.87503112898e-05 -7.77855513392e-05
|
||||
-420 0 0.0373266801818 0.00888730480519 7.93373132824e-05 4.0050086032e-05 4.0050086032e-05 -7.93373132824e-05
|
||||
-415 0 0.0377273343674 0.00909092394396 8.09315481854e-05 4.14074187048e-05 4.14074187048e-05 -8.09315481854e-05
|
||||
-410 0 0.0381360686374 0.00930148015547 8.25695683649e-05 4.28254560779e-05 4.28254560779e-05 -8.25695683649e-05
|
||||
-405 0 0.0385531052638 0.00951928525033 8.42527171775e-05 4.43075480693e-05 4.43075480693e-05 -8.42527171775e-05
|
||||
-400 0 0.0389786733103 0.00974466832758 8.59823675963e-05 4.5857262718e-05 4.5857262718e-05 -8.59823675963e-05
|
||||
-395 0 0.0394130087756 0.00997797690521 8.77599203484e-05 4.74784023799e-05 4.74784023799e-05 -8.77599203484e-05
|
||||
-390 0 0.0398563547266 0.010219578135 8.95868015415e-05 4.91750212543e-05 4.91750212543e-05 -8.95868015415e-05
|
||||
-385 0 0.0403089614182 0.0104698601086 9.14644596945e-05 5.09514443586e-05 5.09514443586e-05 -9.14644596945e-05
|
||||
-380 0 0.0407710863965 0.0107292332622 9.33943620754e-05 5.28122880783e-05 5.28122880783e-05 -9.33943620754e-05
|
||||
-375 0 0.0412429945814 0.0109981318884 9.53779902354e-05 5.47624824318e-05 5.47624824318e-05 -9.53779902354e-05
|
||||
-370 0 0.0417249583257 0.0112770157637 9.74168346095e-05 5.68072952018e-05 5.68072952018e-05 -9.74168346095e-05
|
||||
-365 0 0.0422172574429 0.0115663719022 9.95123880359e-05 5.8952358098e-05 5.8952358098e-05 -9.95123880359e-05
|
||||
-360 0 0.0427201792011 0.0118667164448 0.000101666138022 6.12036951306e-05 6.12036951306e-05 -0.000101666138022
|
||||
-355 0 0.0432340182738 0.0121785966968 0.000103879557563 6.35677533892e-05 6.35677533892e-05 -0.000103879557563
|
||||
-350 0 0.0437590766408 0.0125025933259 0.000106154094277 6.6051436439e-05 6.6051436439e-05 -0.000106154094277
|
||||
-345 0 0.0442956634305 0.0128393227335 0.000108491157607 6.86621405627e-05 6.86621405627e-05 -0.000108491157607
|
||||
-340 0 0.0448440946922 0.0131894396154 0.000110892103772 7.14077940959e-05 7.14077940959e-05 -0.000110892103772
|
||||
-335 0 0.045404693087 0.0135536397275 0.000113358218134 7.42969001219e-05 7.42969001219e-05 -0.000113358218134
|
||||
-330 0 0.0459777874834 0.0139326628737 0.000115890694551 7.73385828147e-05 7.73385828147e-05 -0.000115890694551
|
||||
-325 0 0.046563712441 0.0143272961357 0.000118490611284 8.05426377357e-05 8.05426377357e-05 -0.000118490611284
|
||||
-320 0 0.0471628075641 0.0147383773638 0.000121158902884 8.39195864129e-05 8.39195864129e-05 -0.000121158902884
|
||||
-315 0 0.0477754167033 0.0151667989534 0.000123896327454 8.74807355519e-05 8.74807355519e-05 -0.000123896327454
|
||||
-310 0 0.0484018869805 0.0156135119292 0.000126703428568 9.12382412451e-05 9.12382412451e-05 -0.000126703428568
|
||||
-305 0 0.0490425676083 0.0160795303634 0.00012958049099 9.52051785651e-05 9.52051785651e-05 -0.00012958049099
|
||||
-300 0 0.0496978084706 0.0165659361569 0.000132527489255 9.93956169413e-05 9.93956169413e-05 -0.000132527489255
|
||||
-295 0 0.050367958426 0.0170738842122 0.000135544027977 0.000103824701728 0.000103824701728 -0.000135544027977
|
||||
-290 0 0.051053363289 0.0176046080307 0.00013862927259 0.000108508742378 0.000108508742378 -0.00013862927259
|
||||
-285 0 0.0517543634389 0.018159425768 0.000141781869041 0.000113465307622 0.000113465307622 -0.000141781869041
|
||||
-280 0 0.0524712909946 0.0187397467838 0.00014499985068 0.000118713328042 0.000118713328042 -0.00014499985068
|
||||
-275 0 0.0532044664893 0.0193470787234 0.000148280530362 0.000124273206398 0.000124273206398 -0.000148280530362
|
||||
-270 0 0.0539541949621 0.0199830351711 0.000151620375424 0.000130166935976 0.000130166935976 -0.000151620375424
|
||||
-265 0 0.0547207613783 0.0206493439163 0.000155014862868 0.000136418227182 0.000136418227182 -0.000155014862868
|
||||
-260 0 0.0555044252679 0.0213478558723 0.000158458311629 0.000143052642443 0.000143052642443 -0.000158458311629
|
||||
-255 0 0.0563054144619 0.0220805546909 0.000161943688353 0.000150097739319 0.000150097739319 -0.000161943688353
|
||||
-250 0 0.0571239177824 0.0228495671129 0.000165462382542 0.000157583221469 0.000157583221469 -0.000165462382542
|
||||
-245 0 0.0579600765218 0.0236571740905 0.000169003946287 0.000165541096813 0.000165541096813 -0.000169003946287
|
||||
-240 0 0.0588139745215 0.0245058227173 0.000172555793098 0.00017400584178 0.00017400584178 -0.000172555793098
|
||||
-235 0 0.0596856266289 0.025398138991 0.000176102849501 0.000183014569962 0.000183014569962 -0.000176102849501
|
||||
-230 0 0.0605749652795 0.0263369414259 0.000179627152173 0.000192607202797 0.000192607202797 -0.000179627152173
|
||||
-225 0 0.0614818249121 0.0273252555165 0.000183107382327 0.000202826638885 0.000202826638885 -0.000183107382327
|
||||
-220 0 0.0624059238787 0.0283663290358 0.000186518327906 0.000213718917393 0.000213718917393 -0.000186518327906
|
||||
-215 0 0.0633468434634 0.0294636481225 0.000189830262915 0.000225333369367 0.000225333369367 -0.000189830262915
|
||||
-210 0 0.0643040035665 0.0306209540793 0.000193008231812 0.000237722748859 0.000237722748859 -0.000193008231812
|
||||
-205 0 0.0652766345442 0.0318422607533 0.000196011225492 0.000250943333183 0.000250943333183 -0.000196011225492
|
||||
-200 0 0.0662637446275 0.0331318723138 0.000198791233883 0.00026505497851 0.00026505497851 -0.000198791233883
|
||||
-195 0 0.0672640822612 0.0344944011596 0.000201292158771 0.000280121113008 0.000280121113008 -0.000201292158771
|
||||
-190 0 0.0682760926206 0.0359347855898 0.000203448569174 0.000296208644775 0.000296208644775 -0.000203448569174
|
||||
-185 0 0.0692978674737 0.0374583067425 0.000205184280574 0.000313387755675 0.000313387755675 -0.000205184280574
|
||||
-180 0 0.0703270874585 0.0390706041436 0.000206410738872 0.000331731544615 0.000331731544615 -0.000206410738872
|
||||
-175 0 0.0713609557527 0.0407776890016 0.000207025190316 0.000351315474475 0.000351315474475 -0.000207025190316
|
||||
-170 0 0.0723961220224 0.0425859541308 0.000206908620327 0.000372216565668 0.000372216565668 -0.000206908620327
|
||||
-165 0 0.073428595457 0.0445021790648 0.000205923447788 0.000394512265719 0.000394512265719 -0.000205923447788
|
||||
-160 0 0.0744536456489 0.0465335285306 0.000203910967718 0.00041827890814 0.00041827890814 -0.000203910967718
|
||||
-155 0 0.0754656900607 0.0486875419747 0.000200688545539 0.000443589655023 0.000443589655023 -0.000200688545539
|
||||
-150 0 0.0764581668779 0.0509721112519 0.000196046581738 0.000470511796172 0.000470511796172 -0.000196046581738
|
||||
-145 0 0.0774233921917 0.0533954428908 0.000189745288597 0.000499103253452 0.000499103253452 -0.000189745288597
|
||||
-140 0 0.078352400742 0.05596600053 0.00018151135307 0.000529408113122 0.000529408113122 -0.00018151135307
|
||||
-135 0 0.0792347699266 0.0586924221679 0.0001710346049 0.000561450982651 0.000561450982651 -0.0001710346049
|
||||
-130 0 0.080058427524 0.0615834057877 0.000157964869864 0.000595229944416 0.000595229944416 -0.000157964869864
|
||||
-125 0 0.0808094446677 0.0646475557342 0.000141909268685 0.000630707860821 0.000630707860821 -0.000141909268685
|
||||
-120 0 0.081471817165 0.0678931809708 0.000122430326341 0.000667801780041 0.000667801780041 -0.000122430326341
|
||||
-115 0 0.0820272403893 0.0713280351211 9.90453878431e-05 0.000706370207873 0.000706370207873 -9.90453878431e-05
|
||||
-110 0 0.0824548858487 0.0749589871352 7.12279968253e-05 0.000746198061979 0.000746198061979 -7.12279968253e-05
|
||||
-105 0 0.082731191271 0.0787916107343 3.8412081333e-05 0.00078697922731 0.00078697922731 -3.8412081333e-05
|
||||
-100 0 0.0828296807844 0.0828296807844 4.05284089438e-18 0.000828296807844 0.000828296807844 -4.05284089438e-18
|
||||
-95 0 0.0827208375771 0.0870745658706 -4.46242847431e-05 0.000869601446277 0.000869601446277 4.46242847431e-05
|
||||
-90 0 0.0823720582386 0.091524509154 -9.60754515982e-05 0.000910188488825 0.000910188488825 9.60754515982e-05
|
||||
-85 0 0.0817477255927 0.096173794815 -0.000154938914724 0.000949175333443 0.000949175333443 0.000154938914724
|
||||
-80 0 0.0808094446677 0.101011805835 -0.00022173323232 0.000985481032533 0.000985481032533 0.00022173323232
|
||||
-75 0 0.079516493553 0.106021991404 -0.000296861575931 0.00101781111748 0.00101781111748 0.000296861575931
|
||||
-70 0 0.0778265457035 0.111180779576 -0.000380551661637 0.00104465162018 0.00104465162018 0.000380551661637
|
||||
-65 0 0.0756967205763 0.116456493194 -0.00047278470875 0.00106427726645 0.00106427726645 0.00047278470875
|
||||
-60 0 0.0730850124568 0.121808354095 -0.000573215783975 0.00107477959495 0.00107477959495 0.000573215783975
|
||||
-55 0 0.0699521296452 0.127185690264 -0.000681090356692 0.00107412099263 0.00107412099263 0.000681090356692
|
||||
-50 0 0.0662637446275 0.132527489255 -0.00079516493553 0.00106021991404 0.00106021991404 0.00079516493553
|
||||
-45 0 0.0619931082794 0.137762462843 -0.000913642944843 0.00103107040797 0.00103107040797 0.000913642944843
|
||||
-40 0 0.0571239177824 0.142809794456 -0.00103413989089 0.000984895134178 0.000984895134178 0.00103413989089
|
||||
-35 0 0.0516532530504 0.147580723001 -0.00115369340253 0.000920325221388 0.000920325221388 0.00115369340253
|
||||
-30 0 0.0455943196978 0.151981065659 -0.00126883275 0.000836593021979 0.000836593021979 0.00126883275
|
||||
-25 0 0.0389786733103 0.155914693241 -0.00137571788154 0.000733716203488 0.000733716203488 0.00137571788154
|
||||
-20 0 0.0318575695325 0.159287847662 -0.00147034936304 0.000612645567932 0.000612645567932 0.00147034936304
|
||||
-15 0 0.0243021068316 0.162014045544 -0.00154883843051 0.000475346832892 0.000475346832892 0.00154883843051
|
||||
-10 0 0.016401916987 0.16401916987 -0.00160771265516 0.000324790435386 0.000324790435386 0.00160771265516
|
||||
-5 0 0.0082623122977 0.165246245954 -0.00164422075151 0.000164834160553 0.000164834160553 0.00164422075151
|
||||
0 0 7.52788673411e-16 0.165659361569 -0.00165659361569 -7.66014245849e-18 -7.66014245849e-18 0.00165659361569
|
||||
5 0 -0.0082623122977 0.165246245954 -0.00164422075151 -0.000164834160553 -0.000164834160553 0.00164422075151
|
||||
10 0 -0.016401916987 0.16401916987 -0.00160771265516 -0.000324790435386 -0.000324790435386 0.00160771265516
|
||||
15 0 -0.0243021068316 0.162014045544 -0.00154883843051 -0.000475346832892 -0.000475346832892 0.00154883843051
|
||||
20 0 -0.0318575695325 0.159287847662 -0.00147034936304 -0.000612645567932 -0.000612645567932 0.00147034936304
|
||||
25 0 -0.0389786733103 0.155914693241 -0.00137571788154 -0.000733716203488 -0.000733716203488 0.00137571788154
|
||||
30 0 -0.0455943196978 0.151981065659 -0.00126883275 -0.000836593021979 -0.000836593021979 0.00126883275
|
||||
35 0 -0.0516532530504 0.147580723001 -0.00115369340253 -0.000920325221388 -0.000920325221388 0.00115369340253
|
||||
40 0 -0.0571239177824 0.142809794456 -0.00103413989089 -0.000984895134178 -0.000984895134178 0.00103413989089
|
||||
45 0 -0.0619931082794 0.137762462843 -0.000913642944843 -0.00103107040797 -0.00103107040797 0.000913642944843
|
||||
50 0 -0.0662637446275 0.132527489255 -0.00079516493553 -0.00106021991404 -0.00106021991404 0.00079516493553
|
||||
55 0 -0.0699521296452 0.127185690264 -0.000681090356692 -0.00107412099263 -0.00107412099263 0.000681090356692
|
||||
60 0 -0.0730850124568 0.121808354095 -0.000573215783975 -0.00107477959495 -0.00107477959495 0.000573215783975
|
||||
65 0 -0.0756967205763 0.116456493194 -0.00047278470875 -0.00106427726645 -0.00106427726645 0.00047278470875
|
||||
70 0 -0.0778265457035 0.111180779576 -0.000380551661637 -0.00104465162018 -0.00104465162018 0.000380551661637
|
||||
75 0 -0.079516493553 0.106021991404 -0.000296861575931 -0.00101781111748 -0.00101781111748 0.000296861575931
|
||||
80 0 -0.0808094446677 0.101011805835 -0.00022173323232 -0.000985481032533 -0.000985481032533 0.00022173323232
|
||||
85 0 -0.0817477255927 0.096173794815 -0.000154938914724 -0.000949175333443 -0.000949175333443 0.000154938914724
|
||||
90 0 -0.0823720582386 0.091524509154 -9.60754515982e-05 -0.000910188488825 -0.000910188488825 9.60754515982e-05
|
||||
95 0 -0.0827208375771 0.0870745658706 -4.46242847431e-05 -0.000869601446277 -0.000869601446277 4.46242847431e-05
|
||||
100 0 -0.0828296807844 0.0828296807844 -1.3813386834e-18 -0.000828296807844 -0.000828296807844 1.3813386834e-18
|
||||
105 0 -0.082731191271 0.0787916107343 3.8412081333e-05 -0.00078697922731 -0.00078697922731 -3.8412081333e-05
|
||||
110 0 -0.0824548858487 0.0749589871352 7.12279968253e-05 -0.000746198061979 -0.000746198061979 -7.12279968253e-05
|
||||
115 0 -0.0820272403893 0.0713280351211 9.90453878431e-05 -0.000706370207873 -0.000706370207873 -9.90453878431e-05
|
||||
120 0 -0.081471817165 0.0678931809708 0.000122430326341 -0.000667801780041 -0.000667801780041 -0.000122430326341
|
||||
125 0 -0.0808094446677 0.0646475557342 0.000141909268685 -0.000630707860821 -0.000630707860821 -0.000141909268685
|
||||
130 0 -0.080058427524 0.0615834057877 0.000157964869864 -0.000595229944416 -0.000595229944416 -0.000157964869864
|
||||
135 0 -0.0792347699266 0.0586924221679 0.0001710346049 -0.000561450982651 -0.000561450982651 -0.0001710346049
|
||||
140 0 -0.078352400742 0.05596600053 0.00018151135307 -0.000529408113122 -0.000529408113122 -0.00018151135307
|
||||
145 0 -0.0774233921917 0.0533954428908 0.000189745288597 -0.000499103253452 -0.000499103253452 -0.000189745288597
|
||||
150 0 -0.0764581668779 0.0509721112519 0.000196046581738 -0.000470511796172 -0.000470511796172 -0.000196046581738
|
||||
155 0 -0.0754656900607 0.0486875419747 0.000200688545539 -0.000443589655023 -0.000443589655023 -0.000200688545539
|
||||
160 0 -0.0744536456489 0.0465335285306 0.000203910967718 -0.00041827890814 -0.00041827890814 -0.000203910967718
|
||||
165 0 -0.073428595457 0.0445021790648 0.000205923447788 -0.000394512265719 -0.000394512265719 -0.000205923447788
|
||||
170 0 -0.0723961220224 0.0425859541308 0.000206908620327 -0.000372216565668 -0.000372216565668 -0.000206908620327
|
||||
175 0 -0.0713609557527 0.0407776890016 0.000207025190316 -0.000351315474475 -0.000351315474475 -0.000207025190316
|
||||
180 0 -0.0703270874585 0.0390706041436 0.000206410738872 -0.000331731544615 -0.000331731544615 -0.000206410738872
|
||||
185 0 -0.0692978674737 0.0374583067425 0.000205184280574 -0.000313387755675 -0.000313387755675 -0.000205184280574
|
||||
190 0 -0.0682760926206 0.0359347855898 0.000203448569174 -0.000296208644775 -0.000296208644775 -0.000203448569174
|
||||
195 0 -0.0672640822612 0.0344944011596 0.000201292158771 -0.000280121113008 -0.000280121113008 -0.000201292158771
|
||||
200 0 -0.0662637446275 0.0331318723138 0.000198791233883 -0.00026505497851 -0.00026505497851 -0.000198791233883
|
||||
205 0 -0.0652766345442 0.0318422607533 0.000196011225492 -0.000250943333183 -0.000250943333183 -0.000196011225492
|
||||
210 0 -0.0643040035665 0.0306209540793 0.000193008231812 -0.000237722748859 -0.000237722748859 -0.000193008231812
|
||||
215 0 -0.0633468434634 0.0294636481225 0.000189830262915 -0.000225333369367 -0.000225333369367 -0.000189830262915
|
||||
220 0 -0.0624059238787 0.0283663290358 0.000186518327906 -0.000213718917393 -0.000213718917393 -0.000186518327906
|
||||
225 0 -0.0614818249121 0.0273252555165 0.000183107382327 -0.000202826638885 -0.000202826638885 -0.000183107382327
|
||||
230 0 -0.0605749652795 0.0263369414259 0.000179627152173 -0.000192607202797 -0.000192607202797 -0.000179627152173
|
||||
235 0 -0.0596856266289 0.025398138991 0.000176102849501 -0.000183014569962 -0.000183014569962 -0.000176102849501
|
||||
240 0 -0.0588139745215 0.0245058227173 0.000172555793098 -0.00017400584178 -0.00017400584178 -0.000172555793098
|
||||
245 0 -0.0579600765218 0.0236571740905 0.000169003946287 -0.000165541096813 -0.000165541096813 -0.000169003946287
|
||||
250 0 -0.0571239177824 0.0228495671129 0.000165462382542 -0.000157583221469 -0.000157583221469 -0.000165462382542
|
||||
255 0 -0.0563054144619 0.0220805546909 0.000161943688353 -0.000150097739319 -0.000150097739319 -0.000161943688353
|
||||
260 0 -0.0555044252679 0.0213478558723 0.000158458311629 -0.000143052642443 -0.000143052642443 -0.000158458311629
|
||||
265 0 -0.0547207613783 0.0206493439163 0.000155014862868 -0.000136418227182 -0.000136418227182 -0.000155014862868
|
||||
270 0 -0.0539541949621 0.0199830351711 0.000151620375424 -0.000130166935976 -0.000130166935976 -0.000151620375424
|
||||
275 0 -0.0532044664893 0.0193470787234 0.000148280530362 -0.000124273206398 -0.000124273206398 -0.000148280530362
|
||||
280 0 -0.0524712909946 0.0187397467838 0.00014499985068 -0.000118713328042 -0.000118713328042 -0.00014499985068
|
||||
285 0 -0.0517543634389 0.018159425768 0.000141781869041 -0.000113465307622 -0.000113465307622 -0.000141781869041
|
||||
290 0 -0.051053363289 0.0176046080307 0.00013862927259 -0.000108508742378 -0.000108508742378 -0.00013862927259
|
||||
295 0 -0.050367958426 0.0170738842122 0.000135544027977 -0.000103824701728 -0.000103824701728 -0.000135544027977
|
||||
300 0 -0.0496978084706 0.0165659361569 0.000132527489255 -9.93956169413e-05 -9.93956169413e-05 -0.000132527489255
|
||||
305 0 -0.0490425676083 0.0160795303634 0.00012958049099 -9.52051785651e-05 -9.52051785651e-05 -0.00012958049099
|
||||
310 0 -0.0484018869805 0.0156135119292 0.000126703428568 -9.12382412451e-05 -9.12382412451e-05 -0.000126703428568
|
||||
315 0 -0.0477754167033 0.0151667989534 0.000123896327454 -8.74807355519e-05 -8.74807355519e-05 -0.000123896327454
|
||||
320 0 -0.0471628075641 0.0147383773638 0.000121158902884 -8.3919586413e-05 -8.3919586413e-05 -0.000121158902884
|
||||
325 0 -0.046563712441 0.0143272961357 0.000118490611284 -8.05426377357e-05 -8.05426377357e-05 -0.000118490611284
|
||||
330 0 -0.0459777874834 0.0139326628737 0.000115890694551 -7.73385828147e-05 -7.73385828147e-05 -0.000115890694551
|
||||
335 0 -0.045404693087 0.0135536397275 0.000113358218134 -7.42969001219e-05 -7.42969001219e-05 -0.000113358218134
|
||||
340 0 -0.0448440946922 0.0131894396154 0.000110892103772 -7.14077940958e-05 -7.14077940958e-05 -0.000110892103772
|
||||
345 0 -0.0442956634305 0.0128393227335 0.000108491157607 -6.86621405627e-05 -6.86621405627e-05 -0.000108491157607
|
||||
350 0 -0.0437590766408 0.0125025933259 0.000106154094277 -6.6051436439e-05 -6.6051436439e-05 -0.000106154094277
|
||||
355 0 -0.0432340182738 0.0121785966968 0.000103879557563 -6.35677533892e-05 -6.35677533892e-05 -0.000103879557563
|
||||
360 0 -0.0427201792011 0.0118667164448 0.000101666138022 -6.12036951305e-05 -6.12036951305e-05 -0.000101666138022
|
||||
365 0 -0.0422172574429 0.0115663719022 9.95123880359e-05 -5.8952358098e-05 -5.8952358098e-05 -9.95123880359e-05
|
||||
370 0 -0.0417249583257 0.0112770157637 9.74168346095e-05 -5.68072952018e-05 -5.68072952018e-05 -9.74168346095e-05
|
||||
375 0 -0.0412429945814 0.0109981318884 9.53779902354e-05 -5.47624824318e-05 -5.47624824318e-05 -9.53779902354e-05
|
||||
380 0 -0.0407710863965 0.0107292332622 9.33943620754e-05 -5.28122880783e-05 -5.28122880783e-05 -9.33943620754e-05
|
||||
385 0 -0.0403089614182 0.0104698601086 9.14644596945e-05 -5.09514443586e-05 -5.09514443586e-05 -9.14644596945e-05
|
||||
390 0 -0.0398563547266 0.010219578135 8.95868015415e-05 -4.91750212543e-05 -4.91750212543e-05 -8.95868015415e-05
|
||||
395 0 -0.0394130087756 0.00997797690521 8.77599203484e-05 -4.74784023799e-05 -4.74784023799e-05 -8.77599203484e-05
|
||||
400 0 -0.0389786733103 0.00974466832758 8.59823675963e-05 -4.5857262718e-05 -4.5857262718e-05 -8.59823675963e-05
|
||||
405 0 -0.0385531052638 0.00951928525033 8.42527171775e-05 -4.43075480693e-05 -4.43075480693e-05 -8.42527171775e-05
|
||||
410 0 -0.0381360686374 0.00930148015547 8.25695683649e-05 -4.28254560779e-05 -4.28254560779e-05 -8.25695683649e-05
|
||||
415 0 -0.0377273343674 0.00909092394396 8.09315481854e-05 -4.14074187048e-05 -4.14074187048e-05 -8.09315481854e-05
|
||||
420 0 -0.0373266801818 0.00888730480519 7.93373132824e-05 -4.0050086032e-05 -4.0050086032e-05 -7.93373132824e-05
|
||||
425 0 -0.0369338904481 0.00869032716427 7.77855513392e-05 -3.87503112898e-05 -3.87503112898e-05 -7.77855513392e-05
|
||||
430 0 -0.0365487560157 0.00849971070132 7.62749821273e-05 -3.75051370094e-05 -3.75051370094e-05 -7.62749821273e-05
|
||||
435 0 -0.0361710740532 0.00831518943751 7.48043582351e-05 -3.63117822092e-05 -3.63117822092e-05 -7.48043582351e-05
|
||||
440 0 -0.0358006478832 0.00813651088256 7.33724655225e-05 -3.51676305336e-05 -3.51676305336e-05 -7.33724655225e-05
|
||||
445 0 -0.0354372868156 0.00796343523945 7.19781233457e-05 -3.40702192675e-05 -3.40702192675e-05 -7.19781233457e-05
|
||||
450 0 -0.0350808059793 0.00779573466206 7.06201845857e-05 -3.3017229157e-05 -3.3017229157e-05 -7.06201845857e-05
|
||||
455 0 -0.0347310261554 0.00763319256163 6.92975355123e-05 -3.20064749733e-05 -3.20064749733e-05 -6.92975355123e-05
|
||||
460 0 -0.0343877736109 0.00747560295888 6.80090955104e-05 -3.10358967607e-05 -3.10358967607e-05 -6.80090955104e-05
|
||||
465 0 -0.0340508799335 0.00732276987817 6.67538166924e-05 -3.01035517149e-05 -3.01035517149e-05 -6.67538166924e-05
|
||||
470 0 -0.0337201818698 0.00717450678081 6.55306834159e-05 -2.92076066434e-05 -2.92076066434e-05 -6.55306834159e-05
|
||||
475 0 -0.0333955211651 0.00703063603475 6.43387117238e-05 -2.83463309624e-05 -2.83463309624e-05 -6.43387117238e-05
|
||||
480 0 -0.0330767444064 0.006890988418 6.31769487241e-05 -2.75180901884e-05 -2.75180901884e-05 -6.31769487241e-05
|
||||
485 0 -0.0327637028692 0.00675540265343 6.20444719184e-05 -2.67213398872e-05 -2.67213398872e-05 -6.20444719184e-05
|
||||
490 0 -0.0324562523665 0.00662372497276 6.09403884939e-05 -2.59546200452e-05 -2.59546200452e-05 -6.09403884939e-05
|
||||
495 0 -0.0321542531033 0.00649580870773 5.98638345862e-05 -2.5216549831e-05 -2.5216549831e-05 -5.98638345862e-05
|
||||
500 0 -0.0318575695325 0.00637151390649 5.88139745215e-05 -2.45058227173e-05 -2.45058227173e-05 -5.88139745215e-05
|
||||
4
data/shell2d/config.toml
Normal file
4
data/shell2d/config.toml
Normal file
@ -0,0 +1,4 @@
|
||||
input-file = "data/shell2d/shell2d_model"
|
||||
output-file = "data/shell2d/shell2d_model_out"
|
||||
rho-name = "block"
|
||||
ab-rho = 0.1
|
||||
31
data/shell2d/plot.gnu
Normal file
31
data/shell2d/plot.gnu
Normal file
@ -0,0 +1,31 @@
|
||||
# 设置终端为png格式，输出文件名为"plot.png"
|
||||
set terminal png size 800,600
|
||||
set output "plot.png"
|
||||
|
||||
# 设置数据文件的格式，指定分隔符为逗号
|
||||
set datafile separator ","
|
||||
set key outside
|
||||
|
||||
# 设置子图布局为2行1列
|
||||
set multiplot layout 2,1
|
||||
|
||||
# 第一幅子图
|
||||
set xlabel "angle (deg)"
|
||||
set ylabel "gravity (mGal)"
|
||||
set xrange [10:-10]
|
||||
set key top left inside
|
||||
plot "data/shell2d/shell2d_model_out.csv" using 2:3 with lines title "Vx", \
|
||||
"data/shell2d/shell2d_model_out.csv" using 2:4 with lines title "Vz"
|
||||
|
||||
# 第二幅子图
|
||||
set xlabel "angle (deg)"
|
||||
set ylabel "gradient (Eo)"
|
||||
set xrange [10:-10]
|
||||
set key top left inside
|
||||
plot "data/shell2d/shell2d_model_out.csv" using 2:5 with lines title "Vxx", \
|
||||
"data/shell2d/shell2d_model_out.csv" using 2:6 with lines title "Vxz", \
|
||||
"data/shell2d/shell2d_model_out.csv" using 2:7 with lines title "Vzx", \
|
||||
"data/shell2d/shell2d_model_out.csv" using 2:8 with lines title "Vzz"
|
||||
|
||||
# 结束多图模式
|
||||
unset multiplot
|
||||
49
data/shell2d/shell2d_model.geo
Normal file
49
data/shell2d/shell2d_model.geo
Normal file
@ -0,0 +1,49 @@
|
||||
lc=10000.0;
|
||||
|
||||
deg=10.0;
|
||||
h=300000.0;
|
||||
R=6371000.0;
|
||||
r=R - h;
|
||||
|
||||
deg_bk=1.0;
|
||||
h_bk=75000.0;
|
||||
R_bk=6321000.0;
|
||||
r_bk=R_bk - h_bk;
|
||||
|
||||
Point(1) = {0, 0, 0, lc};
|
||||
Point(2) = {R*Sin(deg*Pi/180.0), R*Cos(deg*Pi/180.0), 0, lc};
|
||||
Point(3) = {r*Sin(deg*Pi/180.0), r*Cos(deg*Pi/180.0), 0, 4*lc};
|
||||
Point(4) = {-1.0*r*Sin(deg*Pi/180.0), r*Cos(deg*Pi/180.0), 0, 4*lc};
|
||||
Point(5) = {-1.0*R*Sin(deg*Pi/180.0), R*Cos(deg*Pi/180.0), 0, lc};
|
||||
|
||||
Point(6) = {R_bk*Sin(deg_bk*Pi/180.0), R_bk*Cos(deg_bk*Pi/180.0), 0, lc};
|
||||
Point(7) = {r_bk*Sin(deg_bk*Pi/180.0), r_bk*Cos(deg_bk*Pi/180.0), 0, 2*lc};
|
||||
Point(8) = {-1.0*r_bk*Sin(deg_bk*Pi/180.0), r_bk*Cos(deg_bk*Pi/180.0), 0, 2*lc};
|
||||
Point(9) = {-1.0*R_bk*Sin(deg_bk*Pi/180.0), R_bk*Cos(deg_bk*Pi/180.0), 0, lc};
|
||||
|
||||
Line(1) = {2, 3};
|
||||
Circle(2) = {3, 1, 4};
|
||||
Line(3) = {4, 5};
|
||||
Circle(4) = {5, 1, 2};
|
||||
|
||||
Physical Curve(1) = {1, 3};
|
||||
Physical Curve(2) = {2};
|
||||
Physical Curve(3) = {4};
|
||||
|
||||
Line(5) = {6, 7};
|
||||
Circle(6) = {7, 1, 8};
|
||||
Line(7) = {8, 9};
|
||||
Circle(8) = {9, 1, 6};
|
||||
|
||||
Curve Loop(1) = {8, 5, 6, 7};
|
||||
Plane Surface(1) = {1};
|
||||
|
||||
Curve Loop(2) = {4, 1, 2, 3};
|
||||
Plane Surface(2) = {1, 2};
|
||||
|
||||
Physical Surface("block", 1) = {1};
|
||||
Physical Surface("shell", 2) = {2};
|
||||
|
||||
Mesh 2;
|
||||
Mesh.MshFileVersion = 2.2;
|
||||
Save "shell2d_model.msh";
|
||||
202
data/shell2d/shell2d_model_out.csv
Normal file
202
data/shell2d/shell2d_model_out.csv
Normal file
@ -0,0 +1,202 @@
|
||||
rad,deg,gx,gz,gxx,gxz,gzx,gzz
|
||||
6.381e+06,10,19.5758,3.46563,1.69621e-05,2.6536e-08,2.6536e-08,-1.6944e-05
|
||||
6.381e+06,9.9,19.7847,3.501,1.73497e-05,9.24959e-08,9.24959e-08,-1.7287e-05
|
||||
6.381e+06,9.8,19.9977,3.53746,1.7749e-05,1.61855e-07,1.61855e-07,-1.76404e-05
|
||||
6.381e+06,9.7,20.2148,3.57505,1.81605e-05,2.34805e-07,2.34805e-07,-1.80045e-05
|
||||
6.381e+06,9.6,20.4363,3.61383,1.85847e-05,3.11549e-07,3.11549e-07,-1.83798e-05
|
||||
6.381e+06,9.5,20.6623,3.65383,1.90221e-05,3.92305e-07,3.92305e-07,-1.87667e-05
|
||||
6.381e+06,9.4,20.8928,3.69513,1.94733e-05,4.77304e-07,4.77304e-07,-1.91657e-05
|
||||
6.381e+06,9.3,21.1281,3.73776,1.99389e-05,5.66793e-07,5.66793e-07,-1.95773e-05
|
||||
6.381e+06,9.2,21.3683,3.78179,2.04194e-05,6.61038e-07,6.61038e-07,-2.00021e-05
|
||||
6.381e+06,9.1,21.6135,3.82729,2.09155e-05,7.6032e-07,7.6032e-07,-2.04406e-05
|
||||
6.381e+06,9,21.8639,3.87431,2.14279e-05,8.64942e-07,8.64942e-07,-2.08934e-05
|
||||
6.381e+06,8.9,22.1197,3.92293,2.19573e-05,9.75227e-07,9.75227e-07,-2.13611e-05
|
||||
6.381e+06,8.8,22.3811,3.97321,2.25044e-05,1.09152e-06,1.09152e-06,-2.18443e-05
|
||||
6.381e+06,8.7,22.6482,4.02524,2.30701e-05,1.2142e-06,1.2142e-06,-2.23439e-05
|
||||
6.381e+06,8.6,22.9212,4.0791,2.36551e-05,1.34365e-06,1.34365e-06,-2.28605e-05
|
||||
6.381e+06,8.5,23.2004,4.13488,2.42604e-05,1.48031e-06,1.48031e-06,-2.33948e-05
|
||||
6.381e+06,8.4,23.4859,4.19265,2.48869e-05,1.62464e-06,1.62464e-06,-2.39478e-05
|
||||
6.381e+06,8.3,23.7781,4.25254,2.55356e-05,1.77712e-06,1.77712e-06,-2.45201e-05
|
||||
6.381e+06,8.2,24.077,4.31463,2.62074e-05,1.9383e-06,1.9383e-06,-2.51129e-05
|
||||
6.381e+06,8.1,24.383,4.37903,2.69036e-05,2.10874e-06,2.10874e-06,-2.5727e-05
|
||||
6.381e+06,8,24.6963,4.44587,2.76254e-05,2.28906e-06,2.28906e-06,-2.63635e-05
|
||||
6.381e+06,7.9,25.0171,4.51527,2.83738e-05,2.47993e-06,2.47993e-06,-2.70234e-05
|
||||
6.381e+06,7.8,25.3459,4.58735,2.91504e-05,2.68206e-06,2.68206e-06,-2.77079e-05
|
||||
6.381e+06,7.7,25.6828,4.66227,2.99564e-05,2.89624e-06,2.89624e-06,-2.84182e-05
|
||||
6.381e+06,7.6,26.0281,4.74017,3.07934e-05,3.12329e-06,3.12329e-06,-2.91556e-05
|
||||
6.381e+06,7.5,26.3823,4.8212,3.16629e-05,3.36414e-06,3.36414e-06,-2.99215e-05
|
||||
6.381e+06,7.4,26.7455,4.90555,3.25667e-05,3.61976e-06,3.61976e-06,-3.07174e-05
|
||||
6.381e+06,7.3,27.1183,4.99339,3.35065e-05,3.89123e-06,3.89123e-06,-3.15447e-05
|
||||
6.381e+06,7.2,27.501,5.08492,3.44842e-05,4.1797e-06,4.1797e-06,-3.24053e-05
|
||||
6.381e+06,7.1,27.8939,5.18035,3.55018e-05,4.48645e-06,4.48645e-06,-3.33007e-05
|
||||
6.381e+06,7,28.2975,5.27991,3.65617e-05,4.81285e-06,4.81285e-06,-3.4233e-05
|
||||
6.381e+06,6.9,28.7122,5.38383,3.76659e-05,5.16039e-06,5.16039e-06,-3.52041e-05
|
||||
6.381e+06,6.8,29.1386,5.49237,3.88171e-05,5.53071e-06,5.53071e-06,-3.62161e-05
|
||||
6.381e+06,6.7,29.577,5.60581,4.00179e-05,5.9256e-06,5.9256e-06,-3.72714e-05
|
||||
6.381e+06,6.6,30.0281,5.72446,4.12711e-05,6.34701e-06,6.34701e-06,-3.83724e-05
|
||||
6.381e+06,6.5,30.4923,5.84863,4.25797e-05,6.79708e-06,6.79708e-06,-3.95217e-05
|
||||
6.381e+06,6.4,30.9703,5.97867,4.3947e-05,7.27815e-06,7.27815e-06,-4.0722e-05
|
||||
6.381e+06,6.3,31.4627,6.11495,4.53764e-05,7.79281e-06,7.79281e-06,-4.19765e-05
|
||||
6.381e+06,6.2,31.9701,6.25789,4.68717e-05,8.34387e-06,8.34387e-06,-4.32883e-05
|
||||
6.381e+06,6.1,32.4932,6.40793,4.8437e-05,8.93446e-06,8.93446e-06,-4.46608e-05
|
||||
6.381e+06,6,33.0328,6.56553,5.00764e-05,9.56802e-06,9.56802e-06,-4.60978e-05
|
||||
6.381e+06,5.9,33.5896,6.73122,5.17947e-05,1.02483e-05,1.02483e-05,-4.76032e-05
|
||||
6.381e+06,5.8,34.1644,6.90556,5.35968e-05,1.09796e-05,1.09796e-05,-4.91813e-05
|
||||
6.381e+06,5.7,34.7582,7.08917,5.54881e-05,1.17665e-05,1.17665e-05,-5.08366e-05
|
||||
6.381e+06,5.6,35.3719,7.28269,5.74743e-05,1.26142e-05,1.26142e-05,-5.25741e-05
|
||||
6.381e+06,5.5,36.0065,7.48688,5.95618e-05,1.35283e-05,1.35283e-05,-5.43991e-05
|
||||
6.381e+06,5.4,36.663,7.70251,6.17572e-05,1.45154e-05,1.45154e-05,-5.63174e-05
|
||||
6.381e+06,5.3,37.3426,7.93046,6.40678e-05,1.55824e-05,1.55824e-05,-5.8335e-05
|
||||
6.381e+06,5.2,38.0465,8.17168,6.65015e-05,1.67375e-05,1.67375e-05,-6.04586e-05
|
||||
6.381e+06,5.1,38.776,8.42722,6.90667e-05,1.79894e-05,1.79894e-05,-6.26954e-05
|
||||
6.381e+06,5,39.5325,8.69824,7.17726e-05,1.93483e-05,1.93483e-05,-6.5053e-05
|
||||
6.381e+06,4.9,40.3174,8.986,7.46291e-05,2.08253e-05,2.08253e-05,-6.75398e-05
|
||||
6.381e+06,4.8,41.1323,9.29191,7.7647e-05,2.24331e-05,2.24331e-05,-7.01647e-05
|
||||
6.381e+06,4.7,41.9789,9.61752,8.08378e-05,2.4186e-05,2.4186e-05,-7.29374e-05
|
||||
6.381e+06,4.6,42.859,9.96454,8.42141e-05,2.61002e-05,2.61002e-05,-7.58682e-05
|
||||
6.381e+06,4.5,43.7746,10.3349,8.77895e-05,2.81941e-05,2.81941e-05,-7.89684e-05
|
||||
6.381e+06,4.4,44.7277,10.7307,9.15785e-05,3.04886e-05,3.04886e-05,-8.22499e-05
|
||||
6.381e+06,4.3,45.7206,11.1543,9.5597e-05,3.30074e-05,3.30074e-05,-8.57254e-05
|
||||
6.381e+06,4.2,46.7557,11.6085,9.98619e-05,3.57777e-05,3.57777e-05,-8.94088e-05
|
||||
6.381e+06,4.1,47.8354,12.096,0.000104391,3.88308e-05,3.88308e-05,-9.33147e-05
|
||||
6.381e+06,4,48.9626,12.6204,0.000109205,4.22023e-05,4.22023e-05,-9.74584e-05
|
||||
6.381e+06,3.9,50.1402,13.1853,0.000114324,4.59335e-05,4.59335e-05,-0.000101856
|
||||
6.381e+06,3.8,51.3713,13.795,0.00011977,5.0072e-05,5.0072e-05,-0.000106525
|
||||
6.381e+06,3.7,52.6592,14.4543,0.000125565,5.46729e-05,5.46729e-05,-0.000111482
|
||||
6.381e+06,3.6,54.0076,15.1686,0.000131735,5.98001e-05,5.98001e-05,-0.000116745
|
||||
6.381e+06,3.5,55.4203,15.9442,0.000138304,6.55284e-05,6.55284e-05,-0.000122332
|
||||
6.381e+06,3.4,56.9012,16.7882,0.000145295,7.19446e-05,7.19446e-05,-0.000128258
|
||||
6.381e+06,3.3,58.4545,17.7086,0.000152733,7.91509e-05,7.91509e-05,-0.000134538
|
||||
6.381e+06,3.2,60.0847,18.7147,0.000160639,8.72667e-05,8.72667e-05,-0.000141184
|
||||
6.381e+06,3.1,61.7963,19.8174,0.00016903,9.64332e-05,9.64332e-05,-0.0001482
|
||||
6.381e+06,3,63.5939,21.0289,0.000177916,0.000106816,0.000106816,-0.000155585
|
||||
6.381e+06,2.9,65.482,22.3636,0.000187298,0.000118613,0.000118613,-0.000163325
|
||||
6.381e+06,2.8,67.465,23.8381,0.00019716,0.000132054,0.000132054,-0.000171388
|
||||
6.381e+06,2.7,69.5468,25.4718,0.000207463,0.000147416,0.000147416,-0.000179716
|
||||
6.381e+06,2.6,71.7307,27.2873,0.000218132,0.000165023,0.000165023,-0.000188219
|
||||
6.381e+06,2.5,74.0186,29.311,0.000229045,0.000185259,0.000185259,-0.000196752
|
||||
6.381e+06,2.4,76.4104,31.5736,0.000240008,0.000208574,0.000208574,-0.000205102
|
||||
6.381e+06,2.3,78.9034,34.1113,0.000250725,0.000235493,0.000235493,-0.000212953
|
||||
6.381e+06,2.2,81.4906,36.9659,0.00026076,0.000266621,0.000266621,-0.000219851
|
||||
6.381e+06,2.1,84.1586,40.1864,0.000269476,0.000302637,0.000302637,-0.000225147
|
||||
6.381e+06,2,86.8858,43.829,0.000275957,0.000344278,0.000344278,-0.000227925
|
||||
6.381e+06,1.9,89.638,47.9579,0.000278906,0.000392295,0.000392295,-0.000226908
|
||||
6.381e+06,1.8,92.364,52.6444,0.000276515,0.000447358,0.000447358,-0.000220335
|
||||
6.381e+06,1.7,94.989,57.9655,0.000266309,0.000509879,0.000509879,-0.00020583
|
||||
6.381e+06,1.6,97.407,63.9991,0.000245006,0.000579704,0.000579704,-0.000180286
|
||||
6.381e+06,1.5,99.4717,70.8151,0.000208463,0.000655604,0.000655604,-0.00013984
|
||||
6.381e+06,1.4,100.989,78.4599,0.000151893,0.00073454,0.00073454,-8.01286e-05
|
||||
6.381e+06,1.3,101.712,86.9319,7.06728e-05,0.000810793,0.000810793,2.88724e-06
|
||||
6.381e+06,1.2,101.358,96.1495,-3.78868e-05,0.000875464,0.000875464,0.000111208
|
||||
6.381e+06,1.1,99.6383,105.921,-0.00017192,0.000917346,0.000917346,0.00024235
|
||||
6.381e+06,1,96.3255,115.935,-0.000322963,0.000926163,0.000926163,0.000387608
|
||||
6.381e+06,0.9,91.3201,125.802,-0.000477057,0.000897199,0.000897199,0.00053342
|
||||
6.381e+06,0.8,84.6836,135.128,-0.000619649,0.000833742,0.000833742,0.000666208
|
||||
6.381e+06,0.7,76.6121,143.595,-0.000740811,0.000744974,0.000744974,0.000777214
|
||||
6.381e+06,0.6,67.3704,150.999,-0.000836994,0.000641429,0.000641429,0.00086386
|
||||
6.381e+06,0.5,57.2301,157.239,-0.000909443,0.000531619,0.000531619,0.000927999
|
||||
6.381e+06,0.4,46.4313,162.291,-0.000961693,0.000420974,0.000420974,0.000973449
|
||||
6.381e+06,0.3,35.1701,166.169,-0.000997711,0.000312213,0.000312213,0.00100425
|
||||
6.381e+06,0.2,23.6007,168.906,-0.00102093,0.000206152,0.000206152,0.0010238
|
||||
6.381e+06,0.1,11.8444,170.533,-0.00103385,0.000102421,0.000102421,0.00103457
|
||||
6.381e+06,0,7.95866e-11,171.073,-0.001038,1.08948e-15,1.08948e-15,0.001038
|
||||
6.381e+06,-0.1,-11.8444,170.533,-0.00103385,-0.000102421,-0.000102421,0.00103457
|
||||
6.381e+06,-0.2,-23.6007,168.906,-0.00102093,-0.000206152,-0.000206152,0.0010238
|
||||
6.381e+06,-0.3,-35.1701,166.169,-0.000997711,-0.000312213,-0.000312213,0.00100425
|
||||
6.381e+06,-0.4,-46.4313,162.291,-0.000961693,-0.000420974,-0.000420974,0.000973449
|
||||
6.381e+06,-0.5,-57.2301,157.239,-0.000909443,-0.000531619,-0.000531619,0.000927999
|
||||
6.381e+06,-0.6,-67.3704,150.999,-0.000836994,-0.000641429,-0.000641429,0.00086386
|
||||
6.381e+06,-0.7,-76.6121,143.595,-0.000740811,-0.000744974,-0.000744974,0.000777214
|
||||
6.381e+06,-0.8,-84.6836,135.128,-0.000619649,-0.000833742,-0.000833742,0.000666208
|
||||
6.381e+06,-0.9,-91.3201,125.802,-0.000477057,-0.000897199,-0.000897199,0.00053342
|
||||
6.381e+06,-1,-96.3255,115.935,-0.000322963,-0.000926163,-0.000926163,0.000387608
|
||||
6.381e+06,-1.1,-99.6383,105.921,-0.00017192,-0.000917346,-0.000917346,0.00024235
|
||||
6.381e+06,-1.2,-101.358,96.1495,-3.78868e-05,-0.000875464,-0.000875464,0.000111208
|
||||
6.381e+06,-1.3,-101.712,86.9319,7.06728e-05,-0.000810793,-0.000810793,2.88724e-06
|
||||
6.381e+06,-1.4,-100.989,78.4599,0.000151893,-0.00073454,-0.00073454,-8.01286e-05
|
||||
6.381e+06,-1.5,-99.4717,70.8151,0.000208463,-0.000655604,-0.000655604,-0.00013984
|
||||
6.381e+06,-1.6,-97.407,63.9991,0.000245006,-0.000579704,-0.000579704,-0.000180286
|
||||
6.381e+06,-1.7,-94.989,57.9655,0.000266309,-0.000509879,-0.000509879,-0.00020583
|
||||
6.381e+06,-1.8,-92.364,52.6444,0.000276515,-0.000447358,-0.000447358,-0.000220335
|
||||
6.381e+06,-1.9,-89.638,47.9579,0.000278906,-0.000392295,-0.000392295,-0.000226908
|
||||
6.381e+06,-2,-86.8858,43.829,0.000275957,-0.000344278,-0.000344278,-0.000227925
|
||||
6.381e+06,-2.1,-84.1586,40.1864,0.000269476,-0.000302637,-0.000302637,-0.000225147
|
||||
6.381e+06,-2.2,-81.4906,36.9659,0.00026076,-0.000266621,-0.000266621,-0.000219851
|
||||
6.381e+06,-2.3,-78.9034,34.1113,0.000250725,-0.000235493,-0.000235493,-0.000212953
|
||||
6.381e+06,-2.4,-76.4104,31.5736,0.000240008,-0.000208574,-0.000208574,-0.000205102
|
||||
6.381e+06,-2.5,-74.0186,29.311,0.000229045,-0.000185259,-0.000185259,-0.000196752
|
||||
6.381e+06,-2.6,-71.7307,27.2873,0.000218132,-0.000165023,-0.000165023,-0.000188219
|
||||
6.381e+06,-2.7,-69.5468,25.4718,0.000207463,-0.000147416,-0.000147416,-0.000179716
|
||||
6.381e+06,-2.8,-67.465,23.8381,0.00019716,-0.000132054,-0.000132054,-0.000171388
|
||||
6.381e+06,-2.9,-65.482,22.3636,0.000187298,-0.000118613,-0.000118613,-0.000163325
|
||||
6.381e+06,-3,-63.5939,21.0289,0.000177916,-0.000106816,-0.000106816,-0.000155585
|
||||
6.381e+06,-3.1,-61.7963,19.8174,0.00016903,-9.64332e-05,-9.64332e-05,-0.0001482
|
||||
6.381e+06,-3.2,-60.0847,18.7147,0.000160639,-8.72667e-05,-8.72667e-05,-0.000141184
|
||||
6.381e+06,-3.3,-58.4545,17.7086,0.000152733,-7.91509e-05,-7.91509e-05,-0.000134538
|
||||
6.381e+06,-3.4,-56.9012,16.7882,0.000145295,-7.19446e-05,-7.19446e-05,-0.000128258
|
||||
6.381e+06,-3.5,-55.4203,15.9442,0.000138304,-6.55284e-05,-6.55284e-05,-0.000122332
|
||||
6.381e+06,-3.6,-54.0076,15.1686,0.000131735,-5.98001e-05,-5.98001e-05,-0.000116745
|
||||
6.381e+06,-3.7,-52.6592,14.4543,0.000125565,-5.46729e-05,-5.46729e-05,-0.000111482
|
||||
6.381e+06,-3.8,-51.3713,13.795,0.00011977,-5.0072e-05,-5.0072e-05,-0.000106525
|
||||
6.381e+06,-3.9,-50.1402,13.1853,0.000114324,-4.59335e-05,-4.59335e-05,-0.000101856
|
||||
6.381e+06,-4,-48.9626,12.6204,0.000109205,-4.22023e-05,-4.22023e-05,-9.74584e-05
|
||||
6.381e+06,-4.1,-47.8354,12.096,0.000104391,-3.88308e-05,-3.88308e-05,-9.33147e-05
|
||||
6.381e+06,-4.2,-46.7557,11.6085,9.98619e-05,-3.57777e-05,-3.57777e-05,-8.94088e-05
|
||||
6.381e+06,-4.3,-45.7206,11.1543,9.5597e-05,-3.30074e-05,-3.30074e-05,-8.57254e-05
|
||||
6.381e+06,-4.4,-44.7277,10.7307,9.15785e-05,-3.04886e-05,-3.04886e-05,-8.22499e-05
|
||||
6.381e+06,-4.5,-43.7746,10.3349,8.77895e-05,-2.81941e-05,-2.81941e-05,-7.89684e-05
|
||||
6.381e+06,-4.6,-42.859,9.96454,8.42141e-05,-2.61002e-05,-2.61002e-05,-7.58682e-05
|
||||
6.381e+06,-4.7,-41.9789,9.61752,8.08378e-05,-2.4186e-05,-2.4186e-05,-7.29374e-05
|
||||
6.381e+06,-4.8,-41.1323,9.29191,7.7647e-05,-2.24331e-05,-2.24331e-05,-7.01647e-05
|
||||
6.381e+06,-4.9,-40.3174,8.986,7.46291e-05,-2.08253e-05,-2.08253e-05,-6.75398e-05
|
||||
6.381e+06,-5,-39.5325,8.69824,7.17726e-05,-1.93483e-05,-1.93483e-05,-6.5053e-05
|
||||
6.381e+06,-5.1,-38.776,8.42722,6.90667e-05,-1.79894e-05,-1.79894e-05,-6.26954e-05
|
||||
6.381e+06,-5.2,-38.0465,8.17168,6.65015e-05,-1.67375e-05,-1.67375e-05,-6.04586e-05
|
||||
6.381e+06,-5.3,-37.3426,7.93046,6.40678e-05,-1.55824e-05,-1.55824e-05,-5.8335e-05
|
||||
6.381e+06,-5.4,-36.663,7.70251,6.17572e-05,-1.45154e-05,-1.45154e-05,-5.63174e-05
|
||||
6.381e+06,-5.5,-36.0065,7.48688,5.95618e-05,-1.35283e-05,-1.35283e-05,-5.43991e-05
|
||||
6.381e+06,-5.6,-35.3719,7.28269,5.74743e-05,-1.26142e-05,-1.26142e-05,-5.25741e-05
|
||||
6.381e+06,-5.7,-34.7582,7.08917,5.54881e-05,-1.17665e-05,-1.17665e-05,-5.08366e-05
|
||||
6.381e+06,-5.8,-34.1644,6.90556,5.35968e-05,-1.09796e-05,-1.09796e-05,-4.91813e-05
|
||||
6.381e+06,-5.9,-33.5896,6.73122,5.17947e-05,-1.02483e-05,-1.02483e-05,-4.76032e-05
|
||||
6.381e+06,-6,-33.0328,6.56553,5.00764e-05,-9.56802e-06,-9.56802e-06,-4.60978e-05
|
||||
6.381e+06,-6.1,-32.4932,6.40793,4.8437e-05,-8.93446e-06,-8.93446e-06,-4.46608e-05
|
||||
6.381e+06,-6.2,-31.9701,6.25789,4.68717e-05,-8.34387e-06,-8.34387e-06,-4.32883e-05
|
||||
6.381e+06,-6.3,-31.4627,6.11495,4.53764e-05,-7.79281e-06,-7.79281e-06,-4.19765e-05
|
||||
6.381e+06,-6.4,-30.9703,5.97867,4.3947e-05,-7.27815e-06,-7.27815e-06,-4.0722e-05
|
||||
6.381e+06,-6.5,-30.4923,5.84863,4.25797e-05,-6.79708e-06,-6.79708e-06,-3.95217e-05
|
||||
6.381e+06,-6.6,-30.0281,5.72446,4.12711e-05,-6.34701e-06,-6.34701e-06,-3.83724e-05
|
||||
6.381e+06,-6.7,-29.577,5.60581,4.00179e-05,-5.9256e-06,-5.9256e-06,-3.72714e-05
|
||||
6.381e+06,-6.8,-29.1386,5.49237,3.88171e-05,-5.53071e-06,-5.53071e-06,-3.62161e-05
|
||||
6.381e+06,-6.9,-28.7122,5.38383,3.76659e-05,-5.16039e-06,-5.16039e-06,-3.52041e-05
|
||||
6.381e+06,-7,-28.2975,5.27991,3.65617e-05,-4.81285e-06,-4.81285e-06,-3.4233e-05
|
||||
6.381e+06,-7.1,-27.8939,5.18035,3.55018e-05,-4.48645e-06,-4.48645e-06,-3.33007e-05
|
||||
6.381e+06,-7.2,-27.501,5.08492,3.44842e-05,-4.1797e-06,-4.1797e-06,-3.24053e-05
|
||||
6.381e+06,-7.3,-27.1183,4.99339,3.35065e-05,-3.89123e-06,-3.89123e-06,-3.15447e-05
|
||||
6.381e+06,-7.4,-26.7455,4.90555,3.25667e-05,-3.61976e-06,-3.61976e-06,-3.07174e-05
|
||||
6.381e+06,-7.5,-26.3823,4.8212,3.16629e-05,-3.36414e-06,-3.36414e-06,-2.99215e-05
|
||||
6.381e+06,-7.6,-26.0281,4.74017,3.07934e-05,-3.12329e-06,-3.12329e-06,-2.91556e-05
|
||||
6.381e+06,-7.7,-25.6828,4.66227,2.99564e-05,-2.89624e-06,-2.89624e-06,-2.84182e-05
|
||||
6.381e+06,-7.8,-25.3459,4.58735,2.91504e-05,-2.68206e-06,-2.68206e-06,-2.77079e-05
|
||||
6.381e+06,-7.9,-25.0171,4.51527,2.83738e-05,-2.47993e-06,-2.47993e-06,-2.70234e-05
|
||||
6.381e+06,-8,-24.6963,4.44587,2.76254e-05,-2.28906e-06,-2.28906e-06,-2.63635e-05
|
||||
6.381e+06,-8.1,-24.383,4.37903,2.69036e-05,-2.10874e-06,-2.10874e-06,-2.5727e-05
|
||||
6.381e+06,-8.2,-24.077,4.31463,2.62074e-05,-1.9383e-06,-1.9383e-06,-2.51129e-05
|
||||
6.381e+06,-8.3,-23.7781,4.25254,2.55356e-05,-1.77712e-06,-1.77712e-06,-2.45201e-05
|
||||
6.381e+06,-8.4,-23.4859,4.19265,2.48869e-05,-1.62464e-06,-1.62464e-06,-2.39478e-05
|
||||
6.381e+06,-8.5,-23.2004,4.13488,2.42604e-05,-1.48031e-06,-1.48031e-06,-2.33948e-05
|
||||
6.381e+06,-8.6,-22.9212,4.0791,2.36551e-05,-1.34365e-06,-1.34365e-06,-2.28605e-05
|
||||
6.381e+06,-8.7,-22.6482,4.02524,2.30701e-05,-1.2142e-06,-1.2142e-06,-2.23439e-05
|
||||
6.381e+06,-8.8,-22.3811,3.97321,2.25044e-05,-1.09152e-06,-1.09152e-06,-2.18443e-05
|
||||
6.381e+06,-8.9,-22.1197,3.92293,2.19573e-05,-9.75227e-07,-9.75227e-07,-2.13611e-05
|
||||
6.381e+06,-9,-21.8639,3.87431,2.14279e-05,-8.64942e-07,-8.64942e-07,-2.08934e-05
|
||||
6.381e+06,-9.1,-21.6135,3.82729,2.09155e-05,-7.6032e-07,-7.6032e-07,-2.04406e-05
|
||||
6.381e+06,-9.2,-21.3683,3.78179,2.04194e-05,-6.61038e-07,-6.61038e-07,-2.00021e-05
|
||||
6.381e+06,-9.3,-21.1281,3.73776,1.99389e-05,-5.66793e-07,-5.66793e-07,-1.95773e-05
|
||||
6.381e+06,-9.4,-20.8928,3.69513,1.94733e-05,-4.77304e-07,-4.77304e-07,-1.91657e-05
|
||||
6.381e+06,-9.5,-20.6623,3.65383,1.90221e-05,-3.92305e-07,-3.92305e-07,-1.87667e-05
|
||||
6.381e+06,-9.6,-20.4363,3.61383,1.85847e-05,-3.11549e-07,-3.11549e-07,-1.83798e-05
|
||||
6.381e+06,-9.7,-20.2148,3.57505,1.81605e-05,-2.34805e-07,-2.34805e-07,-1.80045e-05
|
||||
6.381e+06,-9.8,-19.9977,3.53746,1.7749e-05,-1.61855e-07,-1.61855e-07,-1.76404e-05
|
||||
6.381e+06,-9.9,-19.7847,3.501,1.73497e-05,-9.24959e-08,-9.24959e-08,-1.7287e-05
|
||||
6.381e+06,-10,-19.5758,3.46563,1.69621e-05,-2.6536e-08,-2.6536e-08,-1.6944e-05
|
||||
|
32762
data/stt/mag_obs.csv
Normal file
32762
data/stt/mag_obs.csv
Normal file
File diff suppressed because it is too large
Load Diff
BIN
data/stt/mag_obs.nc
Normal file
BIN
data/stt/mag_obs.nc
Normal file
Binary file not shown.
BIN
data/stt/mag_obs_ren17.nc
Normal file
BIN
data/stt/mag_obs_ren17.nc
Normal file
Binary file not shown.
@ -12,6 +12,8 @@ macro(add_example name switch)
|
||||
endif()
|
||||
endmacro()
|
||||
|
||||
add_example(gobser_tri2d_ex ON)
|
||||
add_example(gobser_tri2d_sph_ex ON)
|
||||
add_example(mobser_dipole_ex OFF)
|
||||
add_example(mobser_block_ex OFF)
|
||||
add_example(mobser_block_gradient_ex OFF)
|
||||
|
||||
79
example/gobser_tri2d_ex.cpp
Normal file
79
example/gobser_tri2d_ex.cpp
Normal file
@ -0,0 +1,79 @@
|
||||
/********************************************************
|
||||
* ██████╗ ██████╗████████╗██╗
|
||||
* ██╔════╝ ██╔════╝╚══██╔══╝██║
|
||||
* ██║ ███╗██║ ██║ ██║
|
||||
* ██║ ██║██║ ██║ ██║
|
||||
* ╚██████╔╝╚██████╗ ██║ ███████╗
|
||||
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
|
||||
* Geophysical Computational Tools & Library (GCTL)
|
||||
*
|
||||
* Copyright (c) 2022 Yi Zhang (yizhang-geo@zju.edu.cn)
|
||||
*
|
||||
* GCTL is distributed under a dual licensing scheme. You can redistribute
|
||||
* it and/or modify it under the terms of the GNU Lesser General Public
|
||||
* License as published by the Free Software Foundation, either version 2
|
||||
* of the License, or (at your option) any later version. You should have
|
||||
* received a copy of the GNU Lesser General Public License along with this
|
||||
* program. If not, see <http://www.gnu.org/licenses/>.
|
||||
*
|
||||
* If the terms and conditions of the LGPL v.2. would prevent you from using
|
||||
* the GCTL, please consider the option to obtain a commercial license for a
|
||||
* fee. These licenses are offered by the GCTL's original author. As a rule,
|
||||
* licenses are provided "as-is", unlimited in time for a one time fee. Please
|
||||
* send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget
|
||||
* to include some description of your company and the realm of its activities.
|
||||
* Also add information on how to contact you by electronic and paper mail.
|
||||
******************************************************/
|
||||
|
||||
#include "gctl/core.h"
|
||||
#include "gctl/io.h"
|
||||
#include "../lib/potential.h"
|
||||
|
||||
using namespace gctl;
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
array<vertex2dc> vert;
|
||||
array<triangle2d> elem;
|
||||
|
||||
gmshio gio;
|
||||
gio.init_file("data/cylinder2d/cylinder2d", Input);
|
||||
gio.set_packed(NotPacked, Input);
|
||||
gio.read_mesh(elem, vert);
|
||||
|
||||
for (size_t i = 0; i < elem.size(); i++)
|
||||
{
|
||||
if (cross(*elem[i].vert[1] - *elem[i].vert[0], *elem[i].vert[2] - *elem[i].vert[0]) < 0.0)
|
||||
{
|
||||
vertex2dc *tmp = elem[i].vert[0]; elem[i].vert[0] = elem[i].vert[1]; elem[i].vert[1] = tmp;
|
||||
}
|
||||
}
|
||||
|
||||
array<double> rho(elem.size(), 1.0);
|
||||
|
||||
array<point2dc> obsp(201);
|
||||
obsp.sequence(point2dc(-500.0, 0.0), point2dc(5.0, 0.0));
|
||||
|
||||
array<double> vx(obsp.size()), vz(obsp.size());
|
||||
array<double> vxx(obsp.size()), vzz(obsp.size());
|
||||
array<double> vxz(obsp.size()), vzx(obsp.size());
|
||||
gobser(vx, elem, obsp, rho, gctl::Vx);
|
||||
gobser(vz, elem, obsp, rho, gctl::Vz);
|
||||
gobser(vxx, elem, obsp, rho, gctl::Txx);
|
||||
gobser(vxz, elem, obsp, rho, gctl::Txz);
|
||||
gobser(vzx, elem, obsp, rho, gctl::Tzx);
|
||||
gobser(vzz, elem, obsp, rho, gctl::Tzz);
|
||||
|
||||
geodsv_io tio;
|
||||
tio.init_table(obsp.size(), 8);
|
||||
tio.set_column_names({"X", "Y", "Vx", "Vz", "Vxx", "Vxz", "Vzx", "Vzz"});
|
||||
tio.fill_column_point2dc(obsp, "X", "Y");
|
||||
tio.fill_column(vx, "Vx", 12);
|
||||
tio.fill_column(vz, "Vz", 12);
|
||||
tio.fill_column(vxx, "Vxx", 12);
|
||||
tio.fill_column(vxz, "Vxz", 12);
|
||||
tio.fill_column(vzx, "Vzx", 12);
|
||||
tio.fill_column(vzz, "Vzz", 12);
|
||||
tio.save_text("data/cylinder2d/gravity_2d");
|
||||
return 0;
|
||||
}
|
||||
114
example/gobser_tri2d_sph_ex.cpp
Normal file
114
example/gobser_tri2d_sph_ex.cpp
Normal file
@ -0,0 +1,114 @@
|
||||
/********************************************************
|
||||
* ██████╗ ██████╗████████╗██╗
|
||||
* ██╔════╝ ██╔════╝╚══██╔══╝██║
|
||||
* ██║ ███╗██║ ██║ ██║
|
||||
* ██║ ██║██║ ██║ ██║
|
||||
* ╚██████╔╝╚██████╗ ██║ ███████╗
|
||||
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
|
||||
* Geophysical Computational Tools & Library (GCTL)
|
||||
*
|
||||
* Copyright (c) 2022 Yi Zhang (yizhang-geo@zju.edu.cn)
|
||||
*
|
||||
* GCTL is distributed under a dual licensing scheme. You can redistribute
|
||||
* it and/or modify it under the terms of the GNU Lesser General Public
|
||||
* License as published by the Free Software Foundation, either version 2
|
||||
* of the License, or (at your option) any later version. You should have
|
||||
* received a copy of the GNU Lesser General Public License along with this
|
||||
* program. If not, see <http://www.gnu.org/licenses/>.
|
||||
*
|
||||
* If the terms and conditions of the LGPL v.2. would prevent you from using
|
||||
* the GCTL, please consider the option to obtain a commercial license for a
|
||||
* fee. These licenses are offered by the GCTL's original author. As a rule,
|
||||
* licenses are provided "as-is", unlimited in time for a one time fee. Please
|
||||
* send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget
|
||||
* to include some description of your company and the realm of its activities.
|
||||
* Also add information on how to contact you by electronic and paper mail.
|
||||
******************************************************/
|
||||
|
||||
#include "gctl/core.h"
|
||||
#include "gctl/geometry.h"
|
||||
#include "gctl/io.h"
|
||||
#include "../lib/potential.h"
|
||||
|
||||
#include "toml.hpp"
|
||||
|
||||
using namespace gctl;
|
||||
int main(int argc, char *argv[]) try
|
||||
{
|
||||
toml::value toml_data = toml::parse(argv[1]);
|
||||
std::string in_file = toml::find<std::string>(toml_data, "input-file");
|
||||
std::string out_file = toml::find<std::string>(toml_data, "output-file");
|
||||
std::string rho_name = toml::find<std::string>(toml_data, "rho-name");
|
||||
double ab_rho = toml::find<double>(toml_data, "ab-rho");
|
||||
|
||||
// read model
|
||||
array<triangle2d> ele;
|
||||
array<vertex2dc> vert;
|
||||
array<gmsh_physical_group> phys;
|
||||
_2i_vector ele_tag;
|
||||
|
||||
gmshio mio;
|
||||
mio.init_file(in_file, Input);
|
||||
mio.set_packed(NotPacked, Input);
|
||||
mio.read_mesh(ele, vert, &ele_tag);
|
||||
mio.read_physical_groups(phys);
|
||||
|
||||
for (size_t i = 0; i < ele.size(); i++)
|
||||
{
|
||||
if (cross(*ele[i].vert[1] - *ele[i].vert[0], *ele[i].vert[2] - *ele[i].vert[0]) < 0.0)
|
||||
{
|
||||
vertex2dc *tmp = ele[i].vert[0]; ele[i].vert[0] = ele[i].vert[1]; ele[i].vert[1] = tmp;
|
||||
}
|
||||
}
|
||||
|
||||
array<double> rho(ele.size());
|
||||
for (size_t i = 0; i < ele.size(); i++)
|
||||
{
|
||||
if (ele_tag[i][0] == mio.physical_name2tag(phys, rho_name)) rho[i] = ab_rho;
|
||||
else rho[i] = 0.0;
|
||||
}
|
||||
|
||||
// calculate gravity
|
||||
array<point2dp> obsp;
|
||||
grid_points_1d(obsp, 100.0, 80.0, -0.1, 6381000.0);
|
||||
|
||||
array<double> gx(obsp.size());
|
||||
array<double> gz(obsp.size());
|
||||
array<double> gxx(obsp.size());
|
||||
array<double> gxz(obsp.size());
|
||||
array<double> gzx(obsp.size());
|
||||
array<double> gzz(obsp.size());
|
||||
gobser(gx, ele, obsp, rho, Vx, FullMsg);
|
||||
gobser(gz, ele, obsp, rho, Vz, FullMsg);
|
||||
gobser(gxx, ele, obsp, rho, Txx, FullMsg);
|
||||
gobser(gxz, ele, obsp, rho, Txz, FullMsg);
|
||||
gobser(gzx, ele, obsp, rho, Tzx, FullMsg);
|
||||
gobser(gzz, ele, obsp, rho, Tzz, FullMsg);
|
||||
|
||||
array<double> r(obsp.size()), d(obsp.size());
|
||||
for (size_t i = 0; i < obsp.size(); i++)
|
||||
{
|
||||
r[i] = obsp[i].rad;
|
||||
d[i] = deg(obsp[i].arc) - 90.0;
|
||||
}
|
||||
|
||||
geodsv_io gio;
|
||||
gio.init_table(obsp.size(), 8);
|
||||
gio.set_column_names({"rad", "deg", "gx", "gz", "gxx", "gxz", "gzx", "gzz"});
|
||||
gio.fill_column(r, "rad");
|
||||
gio.fill_column(d, "deg");
|
||||
gio.fill_column(gx, "gx");
|
||||
gio.fill_column(gz, "gz");
|
||||
gio.fill_column(gxx, "gxx");
|
||||
gio.fill_column(gxz, "gxz");
|
||||
gio.fill_column(gzx, "gzx");
|
||||
gio.fill_column(gzz, "gzz");
|
||||
gio.save_csv(out_file);
|
||||
|
||||
system("gnuplot data/shell2d/plot.gnu");
|
||||
return 0;
|
||||
}
|
||||
catch (std::exception &e)
|
||||
{
|
||||
GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
|
||||
}
|
||||
@ -31,17 +31,36 @@
|
||||
|
||||
using namespace gctl;
|
||||
|
||||
int main(int argc, char const *argv[])
|
||||
bool clockwise(const triangle &t)
|
||||
{
|
||||
if (dot(cross(*t.vert[1] - *t.vert[0], *t.vert[2] - *t.vert[0]), *t.vert[0]) < 0) return true;
|
||||
else return false;
|
||||
}
|
||||
|
||||
int main(int argc, char const *argv[]) try
|
||||
{
|
||||
array<vertex3dc> top_node, btm_node;
|
||||
array<mag_tricone> top_ele, btm_ele;
|
||||
array<magcone_ren17> top_ele, btm_ele;
|
||||
|
||||
gmshio fio;
|
||||
fio.init_file("data/stt/stt_1d_6371200_sph.msh", gctl::Input);
|
||||
fio.set_packed(gctl::Packed, gctl::Input);
|
||||
fio.init_file("data/stt/sph_rect.msh", gctl::Input);
|
||||
fio.set_packed(gctl::NotPacked, gctl::Input);
|
||||
fio.read_mesh(top_ele, top_node);
|
||||
fio.read_mesh(btm_ele, btm_node);
|
||||
|
||||
// 确定三角形顶点排序
|
||||
triangle t;
|
||||
vertex3dc *v = nullptr;
|
||||
for (size_t i = 0; i < top_ele.size(); i++)
|
||||
{
|
||||
t.set(*top_ele[i].vert[0], *top_ele[i].vert[1], *top_ele[i].vert[2]);
|
||||
if (clockwise(t))
|
||||
{
|
||||
v = top_ele[i].vert[0]; top_ele[i].vert[0] = top_ele[i].vert[1]; top_ele[i].vert[1] = v;
|
||||
v = btm_ele[i].vert[0]; btm_ele[i].vert[0] = btm_ele[i].vert[1]; btm_ele[i].vert[1] = v;
|
||||
}
|
||||
}
|
||||
|
||||
vertex3dc origin(point3dc(0.0, 0.0, 0.0), top_node.size());
|
||||
for (size_t i = 0; i < top_ele.size(); i++)
|
||||
{
|
||||
@ -68,49 +87,34 @@ int main(int argc, char const *argv[])
|
||||
mag_B[i] = magkernel_single(md, cen_c);
|
||||
}
|
||||
|
||||
array<magcone_para> top_para, btm_para;
|
||||
// 计算磁体参数
|
||||
array<magcone_para_ren17> top_para, btm_para;
|
||||
array<double> sus(top_ele.size(), 0.01);
|
||||
callink_magnetic_para(top_ele, top_para, mag_B);
|
||||
callink_magnetic_para(btm_ele, btm_para, mag_B);
|
||||
|
||||
// 设置观测点位
|
||||
array<point3ds> obsp(top_ele.size());
|
||||
for (size_t i = 0; i < top_ele.size(); i++)
|
||||
{
|
||||
cen_c = 1.0/3.0*(*top_ele[i].vert[0] + *top_ele[i].vert[1] + *top_ele[i].vert[2]);
|
||||
|
||||
obsp[i] = cen_c.c2s();
|
||||
obsp[i].rad += 200000;
|
||||
}
|
||||
array<point3ds> obsp;
|
||||
grid_points_2d(obsp, 10.0, 55.0, 20.0, 65.0, 0.25, 0.25, 6371200.0 + 100000);
|
||||
|
||||
// 正演计算
|
||||
array<double> obsval(obsp.size());
|
||||
array<point3dc> obsgrad;
|
||||
|
||||
// 正演磁分量数据
|
||||
set_magcone_ren17_tolerance(1e-30);
|
||||
magobser(obsgrad, top_ele, btm_ele, obsp, sus, ShortMsg);
|
||||
|
||||
// 保存网格
|
||||
fio.init_file("data/stt/stt_1d_out.msh", gctl::Output);
|
||||
fio.set_packed(gctl::NotPacked, gctl::Output);
|
||||
fio.save_mesh(top_ele, top_node);
|
||||
geodsv_io fout;
|
||||
fout.init_table(obsp.size(), 6);
|
||||
fout.set_column_names({"rad", "lon", "lat", "Br", "Bt", "Bp"});
|
||||
fout.fill_column_point3ds(obsp, "rad", "lon", "lat");
|
||||
fout.fill_column_point3dc(obsgrad, "Br", "Bt", "Bp", 12);
|
||||
fout.save_csv("data/stt/mag_obs");
|
||||
|
||||
for (int i = 0; i < obsp.size(); ++i)
|
||||
{
|
||||
obsval[i] = obsgrad[i].x;
|
||||
}
|
||||
fio.save_data("Br", obsval, ElemData);
|
||||
|
||||
for (int i = 0; i < obsp.size(); ++i)
|
||||
{
|
||||
obsval[i] = obsgrad[i].y;
|
||||
}
|
||||
fio.save_data("Bt", obsval, ElemData);
|
||||
|
||||
for (int i = 0; i < obsp.size(); ++i)
|
||||
{
|
||||
obsval[i] = obsgrad[i].z;
|
||||
}
|
||||
fio.save_data("Bp", obsval, ElemData);
|
||||
system("xyz2nc -t data/stt/mag_obs.csv -g data/stt/mag_obs_ren17.nc -r 10/55/20/65 -i 0.25/0.25 -c 1,2,3,4,5 -l x,y,Br,Bt,Bp -d ','");
|
||||
return 0;
|
||||
}
|
||||
catch (std::exception &e)
|
||||
{
|
||||
GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
|
||||
}
|
||||
@ -48,6 +48,7 @@
|
||||
#include "potential/mkernel_block.h"
|
||||
#include "potential/mkernel_triangle.h"
|
||||
#include "potential/mkernel_tricone.h"
|
||||
#include "potential/mkernel_tricone_Ren2017.h"
|
||||
#include "potential/mkernel_tetrahedron.h"
|
||||
#include "potential/mkernel_tetrahedron_Ren2017.h"
|
||||
|
||||
|
||||
@ -42,6 +42,12 @@ typedef void (*gkernel_triangle2d_ptr)(gctl::matrix<double> &out_kernel, const g
|
||||
|
||||
void gkernel_triangle2d_vz(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose);
|
||||
void gkernel_triangle2d_vx(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose);
|
||||
void gkernel_triangle2d_vxx(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose);
|
||||
void gkernel_triangle2d_vxz(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose);
|
||||
void gkernel_triangle2d_vzx(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose);
|
||||
void gkernel_triangle2d_vzz(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
@ -65,6 +71,15 @@ void gctl::gkernel(matrix<double> &out_kernel, const array<triangle2d> &ele,
|
||||
case Vz:
|
||||
triangle_kernel = gkernel_triangle2d_vz;
|
||||
break;
|
||||
case Vx:
|
||||
triangle_kernel = gkernel_triangle2d_vx;
|
||||
break;
|
||||
case Txx:
|
||||
triangle_kernel = gkernel_triangle2d_vxx;
|
||||
break;
|
||||
case Txz:
|
||||
triangle_kernel = gkernel_triangle2d_vxz;
|
||||
break;
|
||||
case Tzx:
|
||||
triangle_kernel = gkernel_triangle2d_vzx;
|
||||
break;
|
||||
@ -72,12 +87,59 @@ void gctl::gkernel(matrix<double> &out_kernel, const array<triangle2d> &ele,
|
||||
triangle_kernel = gkernel_triangle2d_vzz;
|
||||
break;
|
||||
default:
|
||||
triangle_kernel = gkernel_triangle2d_vz;
|
||||
throw std::invalid_argument("Invalid gravitational field type!");
|
||||
break;
|
||||
}
|
||||
return triangle_kernel(out_kernel, ele, obsp, verbose);
|
||||
}
|
||||
|
||||
typedef void (*gkernel_triangle2d_ptr2)(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose);
|
||||
|
||||
void gkernel_triangle2d_vz2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose);
|
||||
void gkernel_triangle2d_vx2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose);
|
||||
void gkernel_triangle2d_vxx2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose);
|
||||
void gkernel_triangle2d_vxz2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose);
|
||||
void gkernel_triangle2d_vzx2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose);
|
||||
void gkernel_triangle2d_vzz2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose);
|
||||
|
||||
void gctl::gkernel(matrix<double> &out_kernel, const array<triangle2d> &ele,
|
||||
const array<point2dp> &obsp, gravitational_field_type_e comp_id, verbose_type_e verbose)
|
||||
{
|
||||
gkernel_triangle2d_ptr2 triangle_kernel2;
|
||||
switch (comp_id)
|
||||
{
|
||||
case Vz:
|
||||
triangle_kernel2 = gkernel_triangle2d_vz2;
|
||||
break;
|
||||
case Vx:
|
||||
triangle_kernel2 = gkernel_triangle2d_vx2;
|
||||
break;
|
||||
case Txx:
|
||||
triangle_kernel2 = gkernel_triangle2d_vxx2;
|
||||
break;
|
||||
case Txz:
|
||||
triangle_kernel2 = gkernel_triangle2d_vxz2;
|
||||
break;
|
||||
case Tzx:
|
||||
triangle_kernel2 = gkernel_triangle2d_vzx2;
|
||||
break;
|
||||
case Tzz:
|
||||
triangle_kernel2 = gkernel_triangle2d_vzz2;
|
||||
break;
|
||||
default:
|
||||
throw std::invalid_argument("Invalid gravitational field type!");
|
||||
break;
|
||||
}
|
||||
return triangle_kernel2(out_kernel, ele, obsp, verbose);
|
||||
}
|
||||
|
||||
/**
|
||||
* @brief callback interface of the gravitational observations of 2D triangular elements
|
||||
*
|
||||
@ -92,6 +154,12 @@ typedef void (*gobser_triangle2d_ptr)(gctl::array<double> &out_obs, const gctl::
|
||||
|
||||
void gobser_triangle2d_vz(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
|
||||
void gobser_triangle2d_vx(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
|
||||
void gobser_triangle2d_vxx(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
|
||||
void gobser_triangle2d_vxz(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
|
||||
void gobser_triangle2d_vzx(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
|
||||
void gobser_triangle2d_vzz(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
@ -116,6 +184,15 @@ void gctl::gobser(array<double> &out_obs, const array<triangle2d> &ele, const ar
|
||||
case Vz:
|
||||
triangle_obser = gobser_triangle2d_vz;
|
||||
break;
|
||||
case Vx:
|
||||
triangle_obser = gobser_triangle2d_vx;
|
||||
break;
|
||||
case Txx:
|
||||
triangle_obser = gobser_triangle2d_vxx;
|
||||
break;
|
||||
case Txz:
|
||||
triangle_obser = gobser_triangle2d_vxz;
|
||||
break;
|
||||
case Tzx:
|
||||
triangle_obser = gobser_triangle2d_vzx;
|
||||
break;
|
||||
@ -123,13 +200,63 @@ void gctl::gobser(array<double> &out_obs, const array<triangle2d> &ele, const ar
|
||||
triangle_obser = gobser_triangle2d_vzz;
|
||||
break;
|
||||
default:
|
||||
triangle_obser = gobser_triangle2d_vz;
|
||||
throw std::invalid_argument("Invalid gravitational field type!");
|
||||
break;
|
||||
}
|
||||
return triangle_obser(out_obs, ele, obsp, rho, verbose);
|
||||
}
|
||||
|
||||
typedef void (*gobser_triangle2d_ptr2)(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
|
||||
|
||||
void gobser_triangle2d_vz2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
|
||||
void gobser_triangle2d_vx2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
|
||||
void gobser_triangle2d_vxx2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
|
||||
void gobser_triangle2d_vxz2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
|
||||
void gobser_triangle2d_vzx2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
|
||||
void gobser_triangle2d_vzz2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
|
||||
|
||||
void gctl::gobser(array<double> &out_obs, const array<triangle2d> &ele, const array<point2dp> &obsp,
|
||||
const array<double> &rho, gravitational_field_type_e comp_id, verbose_type_e verbose)
|
||||
{
|
||||
gobser_triangle2d_ptr2 triangle_obser2;
|
||||
switch (comp_id)
|
||||
{
|
||||
case Vz:
|
||||
triangle_obser2 = gobser_triangle2d_vz2;
|
||||
break;
|
||||
case Vx:
|
||||
triangle_obser2 = gobser_triangle2d_vx2;
|
||||
break;
|
||||
case Txx:
|
||||
triangle_obser2 = gobser_triangle2d_vxx2;
|
||||
break;
|
||||
case Txz:
|
||||
triangle_obser2 = gobser_triangle2d_vxz2;
|
||||
break;
|
||||
case Tzx:
|
||||
triangle_obser2 = gobser_triangle2d_vzx2;
|
||||
break;
|
||||
case Tzz:
|
||||
triangle_obser2 = gobser_triangle2d_vzz2;
|
||||
break;
|
||||
default:
|
||||
throw std::invalid_argument("Invalid gravitational field type!");
|
||||
break;
|
||||
}
|
||||
return triangle_obser2(out_obs, ele, obsp, rho, verbose);
|
||||
}
|
||||
|
||||
double gkernel_triangle2d_vz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr);
|
||||
double gkernel_triangle2d_vx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr);
|
||||
double gkernel_triangle2d_vxx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr);
|
||||
double gkernel_triangle2d_vxz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr);
|
||||
double gkernel_triangle2d_vzx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr);
|
||||
double gkernel_triangle2d_vzz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr);
|
||||
|
||||
@ -174,7 +301,79 @@ void gkernel_triangle2d_vz(gctl::matrix<double> &out_kernel, const gctl::array<g
|
||||
#pragma omp parallel for private (j) schedule(guided)
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
out_kernel.at(i, j) = gkernel_triangle2d_vz_sig(ele.get(j), obsp.get(i));
|
||||
out_kernel[i][j] = gkernel_triangle2d_vz_sig(ele.get(j), obsp.get(i));
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
void gkernel_triangle2d_vx(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = ele.size();
|
||||
|
||||
out_kernel.resize(o_size, e_size);
|
||||
|
||||
gctl::progress_bar bar(o_size, "gkernel_vx");
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(i);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);
|
||||
|
||||
#pragma omp parallel for private (j) schedule(guided)
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
out_kernel[i][j] = gkernel_triangle2d_vx_sig(ele.get(j), obsp.get(i));
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
void gkernel_triangle2d_vxx(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = ele.size();
|
||||
|
||||
out_kernel.resize(o_size, e_size);
|
||||
|
||||
gctl::progress_bar bar(o_size, "gkernel_vxx");
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(i);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);
|
||||
|
||||
#pragma omp parallel for private (j) schedule(guided)
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
out_kernel[i][j] = gkernel_triangle2d_vxx_sig(ele.get(j), obsp.get(i));
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
void gkernel_triangle2d_vxz(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = ele.size();
|
||||
|
||||
out_kernel.resize(o_size, e_size);
|
||||
|
||||
gctl::progress_bar bar(o_size, "gkernel_vxz");
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(i);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);
|
||||
|
||||
#pragma omp parallel for private (j) schedule(guided)
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
out_kernel[i][j] = gkernel_triangle2d_vxz_sig(ele.get(j), obsp.get(i));
|
||||
}
|
||||
}
|
||||
return;
|
||||
@ -198,7 +397,7 @@ void gkernel_triangle2d_vzx(gctl::matrix<double> &out_kernel, const gctl::array<
|
||||
#pragma omp parallel for private (j) schedule(guided)
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
out_kernel.at(i, j) = gkernel_triangle2d_vzx_sig(ele.get(j), obsp.get(i));
|
||||
out_kernel[i][j] = gkernel_triangle2d_vzx_sig(ele.get(j), obsp.get(i));
|
||||
}
|
||||
}
|
||||
return;
|
||||
@ -222,12 +421,196 @@ void gkernel_triangle2d_vzz(gctl::matrix<double> &out_kernel, const gctl::array<
|
||||
#pragma omp parallel for private (j) schedule(guided)
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
out_kernel.at(i, j) = gkernel_triangle2d_vzz_sig(ele.get(j), obsp.get(i));
|
||||
out_kernel[i][j] = gkernel_triangle2d_vzz_sig(ele.get(j), obsp.get(i));
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
void gkernel_triangle2d_vz2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = ele.size();
|
||||
|
||||
gctl::point2dc obsc, obsg;
|
||||
out_kernel.resize(o_size, e_size);
|
||||
|
||||
gctl::progress_bar bar(o_size, "gkernel_vz");
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(i);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);
|
||||
|
||||
obsc = obsp[i].p2c();
|
||||
|
||||
#pragma omp parallel for private (j, obsg) schedule(guided)
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
obsg.x = gkernel_triangle2d_vx_sig(ele.get(j), &obsc);
|
||||
obsg.y = gkernel_triangle2d_vz_sig(ele.get(j), &obsc);
|
||||
out_kernel[i][j] = obsg.y*sin(obsp[i].arc) - obsg.x*cos(obsp[i].arc);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void gkernel_triangle2d_vx2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = ele.size();
|
||||
|
||||
gctl::point2dc obsc, obsg;
|
||||
out_kernel.resize(o_size, e_size);
|
||||
|
||||
gctl::progress_bar bar(o_size, "gkernel_vx");
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(i);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);
|
||||
|
||||
obsc = obsp[i].p2c();
|
||||
|
||||
#pragma omp parallel for private (j, obsg) schedule(guided)
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
obsg.x = gkernel_triangle2d_vx_sig(ele.get(j), &obsc);
|
||||
obsg.y = gkernel_triangle2d_vz_sig(ele.get(j), &obsc);
|
||||
out_kernel[i][j] = -1.0*obsg.y*cos(obsp[i].arc) + obsg.x*sin(obsp[i].arc);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void gkernel_triangle2d_vxx2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = ele.size();
|
||||
|
||||
gctl::point2dc obsc;
|
||||
double t[2][2];
|
||||
out_kernel.resize(o_size, e_size);
|
||||
|
||||
gctl::progress_bar bar(o_size, "gkernel_vxx");
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(i);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);
|
||||
|
||||
obsc = obsp[i].p2c();
|
||||
|
||||
#pragma omp parallel for private (j, t) schedule(guided)
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
|
||||
t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
|
||||
t[0][0] = -1.0*t[1][1];
|
||||
t[0][1] = t[1][0];
|
||||
out_kernel[i][j] = sin(obsp[i].arc)*(t[0][0]*sin(obsp[i].arc) - t[1][0]*cos(obsp[i].arc))
|
||||
- cos(obsp[i].arc)*(t[0][1]*sin(obsp[i].arc) - t[1][1]*cos(obsp[i].arc));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void gkernel_triangle2d_vxz2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = ele.size();
|
||||
|
||||
gctl::point2dc obsc;
|
||||
double t[2][2];
|
||||
out_kernel.resize(o_size, e_size);
|
||||
|
||||
gctl::progress_bar bar(o_size, "gkernel_vxz");
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(i);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);
|
||||
|
||||
obsc = obsp[i].p2c();
|
||||
|
||||
#pragma omp parallel for private (j, t) schedule(guided)
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
|
||||
t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
|
||||
t[0][0] = -1.0*t[1][1];
|
||||
t[0][1] = t[1][0];
|
||||
out_kernel[i][j] = -1.0*cos(obsp[i].arc)*(t[0][0]*sin(obsp[i].arc) - t[1][0]*cos(obsp[i].arc))
|
||||
+ sin(obsp[i].arc)*(t[0][1]*sin(obsp[i].arc) - t[1][1]*cos(obsp[i].arc));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void gkernel_triangle2d_vzx2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = ele.size();
|
||||
|
||||
gctl::point2dc obsc;
|
||||
double t[2][2];
|
||||
out_kernel.resize(o_size, e_size);
|
||||
|
||||
gctl::progress_bar bar(o_size, "gkernel_vzx");
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(i);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);
|
||||
|
||||
obsc = obsp[i].p2c();
|
||||
|
||||
#pragma omp parallel for private (j, t) schedule(guided)
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
|
||||
t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
|
||||
t[0][0] = -1.0*t[1][1];
|
||||
t[0][1] = t[1][0];
|
||||
out_kernel[i][j] = sin(obsp[i].arc)*(t[1][0]*sin(obsp[i].arc) - t[0][0]*cos(obsp[i].arc))
|
||||
- cos(obsp[i].arc)*(t[1][1]*sin(obsp[i].arc) - t[0][1]*cos(obsp[i].arc));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void gkernel_triangle2d_vzz2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = ele.size();
|
||||
|
||||
gctl::point2dc obsc;
|
||||
double t[2][2];
|
||||
out_kernel.resize(o_size, e_size);
|
||||
|
||||
gctl::progress_bar bar(o_size, "gkernel_vzz");
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(i);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);
|
||||
|
||||
obsc = obsp[i].p2c();
|
||||
|
||||
#pragma omp parallel for private (j, t) schedule(guided)
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
|
||||
t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
|
||||
t[0][0] = -1.0*t[1][1];
|
||||
t[0][1] = t[1][0];
|
||||
out_kernel[i][j] = -1.0*cos(obsp[i].arc)*(t[1][0]*sin(obsp[i].arc) - t[0][0]*cos(obsp[i].arc))
|
||||
+ sin(obsp[i].arc)*(t[1][1]*sin(obsp[i].arc) - t[0][1]*cos(obsp[i].arc));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void gobser_triangle2d_vz(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
|
||||
{
|
||||
@ -246,7 +629,79 @@ void gobser_triangle2d_vz(gctl::array<double> &out_obs, const gctl::array<gctl::
|
||||
#pragma omp parallel for private (i) schedule(guided)
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
out_obs.at(i) += gkernel_triangle2d_vz_sig(ele.get(j), obsp.get(i)) * rho[j];
|
||||
out_obs[i] += gkernel_triangle2d_vz_sig(ele.get(j), obsp.get(i)) * rho[j];
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
void gobser_triangle2d_vx(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = ele.size();
|
||||
|
||||
out_obs.resize(o_size, 0.0);
|
||||
|
||||
gctl::progress_bar bar(e_size, "gobser_vx");
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(j);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
|
||||
|
||||
#pragma omp parallel for private (i) schedule(guided)
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
out_obs[i] += gkernel_triangle2d_vx_sig(ele.get(j), obsp.get(i)) * rho[j];
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
void gobser_triangle2d_vxx(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = ele.size();
|
||||
|
||||
out_obs.resize(o_size, 0.0);
|
||||
|
||||
gctl::progress_bar bar(e_size, "gobser_vxx");
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(j);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
|
||||
|
||||
#pragma omp parallel for private (i) schedule(guided)
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
out_obs[i] += gkernel_triangle2d_vxx_sig(ele.get(j), obsp.get(i)) * rho[j];
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
void gobser_triangle2d_vxz(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = ele.size();
|
||||
|
||||
out_obs.resize(o_size, 0.0);
|
||||
|
||||
gctl::progress_bar bar(e_size, "gobser_vxz");
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(j);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
|
||||
|
||||
#pragma omp parallel for private (i) schedule(guided)
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
out_obs[i] += gkernel_triangle2d_vxz_sig(ele.get(j), obsp.get(i)) * rho[j];
|
||||
}
|
||||
}
|
||||
return;
|
||||
@ -270,7 +725,7 @@ void gobser_triangle2d_vzx(gctl::array<double> &out_obs, const gctl::array<gctl:
|
||||
#pragma omp parallel for private (i) schedule(guided)
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
out_obs.at(i) += gkernel_triangle2d_vzx_sig(ele.get(j), obsp.get(i)) * rho[j];
|
||||
out_obs[i] += gkernel_triangle2d_vzx_sig(ele.get(j), obsp.get(i)) * rho[j];
|
||||
}
|
||||
}
|
||||
return;
|
||||
@ -294,12 +749,192 @@ void gobser_triangle2d_vzz(gctl::array<double> &out_obs, const gctl::array<gctl:
|
||||
#pragma omp parallel for private (i) schedule(guided)
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
out_obs.at(i) += gkernel_triangle2d_vzz_sig(ele.get(j), obsp.get(i)) * rho[j];
|
||||
out_obs[i] += gkernel_triangle2d_vzz_sig(ele.get(j), obsp.get(i)) * rho[j];
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
void gobser_triangle2d_vz2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = ele.size();
|
||||
|
||||
gctl::point2dc obsc, obsg;
|
||||
out_obs.resize(o_size, 0.0);
|
||||
|
||||
gctl::progress_bar bar(e_size, "gobser_vz");
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(j);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
|
||||
|
||||
#pragma omp parallel for private (i, obsg, obsc) schedule(guided)
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
obsc = obsp[i].p2c();
|
||||
obsg.x = gkernel_triangle2d_vx_sig(ele.get(j), &obsc); // 沿x轴正方向的重力值
|
||||
obsg.y = gkernel_triangle2d_vz_sig(ele.get(j), &obsc); // 沿z轴负方向的重力值
|
||||
out_obs[i] += (obsg.y*sin(obsp[i].arc) - obsg.x*cos(obsp[i].arc))*rho[j]; // 沿半径负方向的重力值
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
void gobser_triangle2d_vx2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = ele.size();
|
||||
|
||||
gctl::point2dc obsc, obsg;
|
||||
out_obs.resize(o_size, 0.0);
|
||||
|
||||
gctl::progress_bar bar(e_size, "gobser_vx");
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(j);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
|
||||
|
||||
#pragma omp parallel for private (i, obsg, obsc) schedule(guided)
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
obsc = obsp[i].p2c();
|
||||
obsg.x = gkernel_triangle2d_vx_sig(ele.get(j), &obsc); // 沿x轴正方向的重力值
|
||||
obsg.y = gkernel_triangle2d_vz_sig(ele.get(j), &obsc); // 沿z轴负方向的重力值
|
||||
out_obs[i] += (-1.0*obsg.y*cos(obsp[i].arc) + obsg.x*sin(obsp[i].arc))*rho[j]; // 沿弧度负方向的重力值
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
void gobser_triangle2d_vxx2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = ele.size();
|
||||
|
||||
gctl::point2dc obsc;
|
||||
double t[2][2];
|
||||
out_obs.resize(o_size, 0.0);
|
||||
|
||||
gctl::progress_bar bar(e_size, "gobser_vxx");
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(j);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
|
||||
|
||||
#pragma omp parallel for private (i, t, obsc) schedule(guided)
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
obsc = obsp[i].p2c();
|
||||
t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
|
||||
t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
|
||||
t[0][0] = -1.0*t[1][1];
|
||||
t[0][1] = t[1][0];
|
||||
out_obs[i] += (sin(obsp[i].arc)*(t[0][0]*sin(obsp[i].arc) - t[1][0]*cos(obsp[i].arc))
|
||||
- cos(obsp[i].arc)*(t[0][1]*sin(obsp[i].arc) - t[1][1]*cos(obsp[i].arc)))*rho[j];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void gobser_triangle2d_vxz2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = ele.size();
|
||||
|
||||
gctl::point2dc obsc;
|
||||
double t[2][2];
|
||||
out_obs.resize(o_size, 0.0);
|
||||
|
||||
gctl::progress_bar bar(e_size, "gobser_vxz");
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(j);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
|
||||
|
||||
#pragma omp parallel for private (i, t, obsc) schedule(guided)
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
obsc = obsp[i].p2c();
|
||||
t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
|
||||
t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
|
||||
t[0][0] = -1.0*t[1][1];
|
||||
t[0][1] = t[1][0];
|
||||
out_obs[i] += (-1.0*cos(obsp[i].arc)*(t[0][0]*sin(obsp[i].arc) - t[1][0]*cos(obsp[i].arc))
|
||||
+ sin(obsp[i].arc)*(t[0][1]*sin(obsp[i].arc) - t[1][1]*cos(obsp[i].arc)))*rho[j];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void gobser_triangle2d_vzx2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = ele.size();
|
||||
|
||||
gctl::point2dc obsc;
|
||||
double t[2][2];
|
||||
out_obs.resize(o_size, 0.0);
|
||||
|
||||
gctl::progress_bar bar(e_size, "gobser_vzx");
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(j);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
|
||||
|
||||
#pragma omp parallel for private (i, t, obsc) schedule(guided)
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
obsc = obsp[i].p2c();
|
||||
t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
|
||||
t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
|
||||
t[0][0] = -1.0*t[1][1];
|
||||
t[0][1] = t[1][0];
|
||||
out_obs[i] += (sin(obsp[i].arc)*(t[1][0]*sin(obsp[i].arc) - t[0][0]*cos(obsp[i].arc))
|
||||
- cos(obsp[i].arc)*(t[1][1]*sin(obsp[i].arc) - t[0][1]*cos(obsp[i].arc)))*rho[j];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void gobser_triangle2d_vzz2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele,
|
||||
const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = ele.size();
|
||||
|
||||
gctl::point2dc obsc;
|
||||
double t[2][2];
|
||||
out_obs.resize(o_size, 0.0);
|
||||
|
||||
gctl::progress_bar bar(e_size, "gobser_vzz");
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(j);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
|
||||
|
||||
#pragma omp parallel for private (i, t, obsc) schedule(guided)
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
obsc = obsp[i].p2c();
|
||||
t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
|
||||
t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
|
||||
t[0][0] = -1.0*t[1][1];
|
||||
t[0][1] = t[1][0];
|
||||
out_obs[i] += (-1.0*cos(obsp[i].arc)*(t[1][0]*sin(obsp[i].arc) - t[0][0]*cos(obsp[i].arc))
|
||||
+ sin(obsp[i].arc)*(t[1][1]*sin(obsp[i].arc) - t[0][1]*cos(obsp[i].arc)))*rho[j];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
double gkernel_triangle2d_vz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr)
|
||||
{
|
||||
double sum;
|
||||
@ -317,17 +952,59 @@ double gkernel_triangle2d_vz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_p
|
||||
|
||||
Ba = atan2(ele_ptr->vert[n]->y - op_ptr->y, ele_ptr->vert[n]->x - op_ptr->x);
|
||||
Bb = atan2(ele_ptr->vert[(n+1)%3]->y - op_ptr->y, ele_ptr->vert[(n+1)%3]->x - op_ptr->x);
|
||||
if (ele_ptr->vert[n]->y - op_ptr->y > 0.0 && ele_ptr->vert[n]->x - op_ptr->x < 0.0) Ba -= 2*GCTL_Pi;
|
||||
if (ele_ptr->vert[(n+1)%3]->y - op_ptr->y > 0.0 && ele_ptr->vert[(n+1)%3]->x - op_ptr->x < 0.0) Bb -= 2*GCTL_Pi;
|
||||
B = (ele_ptr->vert[(n+1)%3]->x - ele_ptr->vert[n]->x)*(Ba - Bb);
|
||||
|
||||
Ca = pow(ele_ptr->vert[(n+1)%3]->x - op_ptr->x,2) + pow(ele_ptr->vert[(n+1)%3]->y - op_ptr->y,2);
|
||||
Cb = pow(ele_ptr->vert[n]->x - op_ptr->x,2) + pow(ele_ptr->vert[n]->y - op_ptr->y,2);
|
||||
C = 0.5*(ele_ptr->vert[(n+1)%3]->y - ele_ptr->vert[n]->y)*(log(Ca) - log(Cb));
|
||||
|
||||
sum += A*(B+C);
|
||||
sum += A*(B + C);
|
||||
}
|
||||
return -2.0e+8*GCTL_G0*sum;
|
||||
}
|
||||
|
||||
double gkernel_triangle2d_vx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr)
|
||||
{
|
||||
double sum;
|
||||
double A, B, C;
|
||||
double Aa, Ab, Ba, Bb, Ca, Cb;
|
||||
|
||||
sum = 0.0;
|
||||
for (int n = 0; n < 3; n++)
|
||||
{
|
||||
Aa = (ele_ptr->vert[n]->y - op_ptr->y)*(ele_ptr->vert[(n+1)%3]->x - op_ptr->x) -
|
||||
(ele_ptr->vert[(n+1)%3]->y - op_ptr->y)*(ele_ptr->vert[n]->x - op_ptr->x);
|
||||
Ab = pow(ele_ptr->vert[(n+1)%3]->y - ele_ptr->vert[n]->y,2) +
|
||||
pow(ele_ptr->vert[(n+1)%3]->x - ele_ptr->vert[n]->x,2);
|
||||
A = Aa/Ab;
|
||||
|
||||
Ba = atan2(ele_ptr->vert[n]->y - op_ptr->y, ele_ptr->vert[n]->x - op_ptr->x);
|
||||
Bb = atan2(ele_ptr->vert[(n+1)%3]->y - op_ptr->y, ele_ptr->vert[(n+1)%3]->x - op_ptr->x);
|
||||
if (ele_ptr->vert[n]->y - op_ptr->y > 0.0 && ele_ptr->vert[n]->x - op_ptr->x < 0.0) Ba -= 2*GCTL_Pi;
|
||||
if (ele_ptr->vert[(n+1)%3]->y - op_ptr->y > 0.0 && ele_ptr->vert[(n+1)%3]->x - op_ptr->x < 0.0) Bb -= 2*GCTL_Pi;
|
||||
B = (ele_ptr->vert[(n+1)%3]->y - ele_ptr->vert[n]->y)*(Ba - Bb);
|
||||
|
||||
Ca = pow(ele_ptr->vert[(n+1)%3]->x - op_ptr->x,2) + pow(ele_ptr->vert[(n+1)%3]->y - op_ptr->y,2);
|
||||
Cb = pow(ele_ptr->vert[n]->x - op_ptr->x,2) + pow(ele_ptr->vert[n]->y - op_ptr->y,2);
|
||||
C = 0.5*(ele_ptr->vert[(n+1)%3]->x - ele_ptr->vert[n]->x)*(log(Cb) - log(Ca));
|
||||
|
||||
sum += A*(B + C);
|
||||
}
|
||||
return 2.0e+8*GCTL_G0*sum;
|
||||
}
|
||||
|
||||
double gkernel_triangle2d_vxx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr)
|
||||
{
|
||||
return -1.0*gkernel_triangle2d_vzz_sig(ele_ptr, op_ptr);
|
||||
}
|
||||
|
||||
double gkernel_triangle2d_vxz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr)
|
||||
{
|
||||
return gkernel_triangle2d_vzx_sig(ele_ptr, op_ptr);
|
||||
}
|
||||
|
||||
double gkernel_triangle2d_vzx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr)
|
||||
{
|
||||
double sum;
|
||||
@ -348,10 +1025,12 @@ double gkernel_triangle2d_vzx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_
|
||||
Ax= Aa_x/Ab;
|
||||
|
||||
Ba = atan2(ele_ptr->vert[n]->y - op_ptr->y, ele_ptr->vert[n]->x - op_ptr->x);
|
||||
if (ele_ptr->vert[n]->y - op_ptr->y > 0.0 && ele_ptr->vert[n]->x - op_ptr->x < 0.0) Ba -= 2*GCTL_Pi;
|
||||
Ba_x = (ele_ptr->vert[n]->y - op_ptr->y)/((ele_ptr->vert[n]->x - op_ptr->x)*(ele_ptr->vert[n]->x - op_ptr->x)
|
||||
+ (ele_ptr->vert[n]->y - op_ptr->y)*(ele_ptr->vert[n]->y - op_ptr->y));
|
||||
|
||||
Bb = atan2(ele_ptr->vert[(n+1)%3]->y - op_ptr->y, ele_ptr->vert[(n+1)%3]->x - op_ptr->x);
|
||||
if (ele_ptr->vert[(n+1)%3]->y - op_ptr->y > 0.0 && ele_ptr->vert[(n+1)%3]->x - op_ptr->x < 0.0) Bb -= 2*GCTL_Pi;
|
||||
Bb_x = (ele_ptr->vert[(n+1)%3]->y - op_ptr->y)/((ele_ptr->vert[(n+1)%3]->x - op_ptr->x)*(ele_ptr->vert[(n+1)%3]->x - op_ptr->x)
|
||||
+ (ele_ptr->vert[(n+1)%3]->y - op_ptr->y)*(ele_ptr->vert[(n+1)%3]->y - op_ptr->y));
|
||||
|
||||
@ -391,10 +1070,12 @@ double gkernel_triangle2d_vzz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_
|
||||
Ay= Aa_y/Ab;
|
||||
|
||||
Ba = atan2(ele_ptr->vert[n]->y - op_ptr->y, ele_ptr->vert[n]->x - op_ptr->x);
|
||||
if (ele_ptr->vert[n]->y - op_ptr->y > 0.0 && ele_ptr->vert[n]->x - op_ptr->x < 0.0) Ba -= 2*GCTL_Pi;
|
||||
Ba_y = -1.0*(ele_ptr->vert[n]->x - op_ptr->x)/((ele_ptr->vert[n]->x - op_ptr->x)*(ele_ptr->vert[n]->x - op_ptr->x)
|
||||
+ (ele_ptr->vert[n]->y - op_ptr->y)*(ele_ptr->vert[n]->y - op_ptr->y));
|
||||
|
||||
Bb = atan2(ele_ptr->vert[(n+1)%3]->y - op_ptr->y, ele_ptr->vert[(n+1)%3]->x - op_ptr->x);
|
||||
if (ele_ptr->vert[(n+1)%3]->y - op_ptr->y > 0.0 && ele_ptr->vert[(n+1)%3]->x - op_ptr->x < 0.0) Bb -= 2*GCTL_Pi;
|
||||
Bb_y = -1.0*(ele_ptr->vert[(n+1)%3]->x - op_ptr->x)/((ele_ptr->vert[(n+1)%3]->x - op_ptr->x)*(ele_ptr->vert[(n+1)%3]->x - op_ptr->x)
|
||||
+ (ele_ptr->vert[(n+1)%3]->y - op_ptr->y)*(ele_ptr->vert[(n+1)%3]->y - op_ptr->y));
|
||||
|
||||
|
||||
@ -56,6 +56,31 @@ namespace gctl
|
||||
*/
|
||||
void gobser(array<double> &out_obs, const array<triangle2d> &ele, const array<point2dc> &obsp,
|
||||
const array<double> &rho, gravitational_field_type_e comp_id = Vz, verbose_type_e verbose = FullMsg);
|
||||
|
||||
/**
|
||||
* @brief 计算极坐标下的二维三角形单元体正演核矩阵
|
||||
*
|
||||
* @param out_kernel 输出的核矩阵
|
||||
* @param[in] ele 二维三角形单元体数组
|
||||
* @param[in] obsp 二维极坐标下的观测点数组
|
||||
* @param[in] comp_id 待计算的重力场分量 默认为重力值
|
||||
* @param[in] verbose 返回信息层级
|
||||
*/
|
||||
void gkernel(matrix<double> &out_kernel, const array<triangle2d> &ele, const array<point2dp> &obsp,
|
||||
gravitational_field_type_e comp_id = Vz, verbose_type_e verbose = FullMsg);
|
||||
|
||||
/**
|
||||
* @brief 计算极坐标下的二维三角形单元体重力观测值正演
|
||||
*
|
||||
* @param out_obs 输出的重力观测值
|
||||
* @param[in] ele 二维三角形单元体数组
|
||||
* @param[in] obsp 二维极坐标下的观测点数组
|
||||
* @param[in] rho 三角形单元体密度数组
|
||||
* @param[in] comp_id 待计算的重力场分量 默认为重力值
|
||||
* @param[in] verbose 返回信息层级
|
||||
*/
|
||||
void gobser(array<double> &out_obs, const array<triangle2d> &ele, const array<point2dp> &obsp,
|
||||
const array<double> &rho, gravitational_field_type_e comp_id = Vz, verbose_type_e verbose = FullMsg);
|
||||
}
|
||||
|
||||
#endif // _GCTL_GRAV_KERNEL_TRIANGLE2D_H
|
||||
@ -133,7 +133,7 @@ void gctl::read_IGRF_table(std::string file, array<IGRF_para> &IGRFs, int head_r
|
||||
}
|
||||
infile.close();
|
||||
|
||||
IGRFs.import_vector(vec_para);
|
||||
IGRFs.input(vec_para);
|
||||
destroy_vector(vec_para);
|
||||
return;
|
||||
}
|
||||
|
||||
@ -31,6 +31,7 @@
|
||||
#include "gctl/core.h"
|
||||
#include "gctl/utility.h"
|
||||
#include "gctl/geometry.h"
|
||||
#include "gctl/maths.h"
|
||||
#include "gctl/algorithm.h"
|
||||
#include "gctl_potential_config.h"
|
||||
|
||||
|
||||
@ -44,9 +44,9 @@ gctl::gm_regular_grid::gm_regular_grid(std::string in_name, std::string in_info,
|
||||
|
||||
void gctl::gm_regular_grid::gradient(std::string datname, std::string gradname, gravitational_field_type_e c_type, int order)
|
||||
{
|
||||
meshdata *data_ptr = get_data(datname);
|
||||
mesh_data_value_e value_type = data_ptr->get_valtype();
|
||||
mesh_data_type_e data_type = data_ptr->get_dattype();
|
||||
meshdata &data = get_data(datname);
|
||||
mesh_data_value_e value_type = data.valtype_;
|
||||
mesh_data_type_e data_type = data.loctype_;
|
||||
|
||||
if(value_type != Scalar)
|
||||
{
|
||||
@ -59,16 +59,13 @@ void gctl::gm_regular_grid::gradient(std::string datname, std::string gradname,
|
||||
}
|
||||
|
||||
// 检查是否存在与梯度数据同名的数据 若有则检查是否需要重新建立数据
|
||||
meshdata *new_data_ptr = add_data(gradname, data_type, true, value_type);
|
||||
|
||||
array<double> *inval_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
array<double> *outval_ptr = (array<double>*) new_data_ptr->get_datval_ptr();
|
||||
meshdata &new_data = add_data(data_type, value_type, gradname, 0.0);
|
||||
gctl::array<double> ingrid_ex;
|
||||
|
||||
int M = -1, N = -1;
|
||||
int ori_row, ori_col;
|
||||
if (data_type == NodeData) gctl::cosine_extrapolate_2d(*inval_ptr, ingrid_ex, rg_ynum, rg_xnum, M, N, ori_row, ori_col);
|
||||
else gctl::cosine_extrapolate_2d(*inval_ptr, ingrid_ex, rg_ynum-1, rg_xnum-1, M, N, ori_row, ori_col);
|
||||
if (data_type == NodeData) gctl::cosine_extrapolate_2d(data.datval_, ingrid_ex, rg_ynum, rg_xnum, M, N, ori_row, ori_col);
|
||||
else gctl::cosine_extrapolate_2d(data.datval_, ingrid_ex, rg_ynum-1, rg_xnum-1, M, N, ori_row, ori_col);
|
||||
|
||||
fftw_complex *in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*M*N);
|
||||
fftw_complex *out= (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*M*N);
|
||||
@ -173,7 +170,7 @@ void gctl::gm_regular_grid::gradient(std::string datname, std::string gradname,
|
||||
{
|
||||
for (int j = 0; j < rg_xnum; j++)
|
||||
{
|
||||
outval_ptr->at(i*rg_xnum+j) = ingrid_ex[(i+ori_row)*N+j+ori_col];
|
||||
new_data.datval_[i*rg_xnum+j] = ingrid_ex[(i+ori_row)*N+j+ori_col];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -183,7 +180,7 @@ void gctl::gm_regular_grid::gradient(std::string datname, std::string gradname,
|
||||
{
|
||||
for (int j = 0; j < rg_xnum-1; j++)
|
||||
{
|
||||
outval_ptr->at(i*(rg_xnum-1)+j) = ingrid_ex[(i+ori_row)*N+j+ori_col];
|
||||
new_data.datval_[i*(rg_xnum-1)+j] = ingrid_ex[(i+ori_row)*N+j+ori_col];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -196,9 +193,9 @@ void gctl::gm_regular_grid::gradient(std::string datname, std::string gradname,
|
||||
|
||||
void gctl::gm_regular_grid::rtp(std::string datname, std::string rtpname, double inc, double dec)
|
||||
{
|
||||
meshdata *data_ptr = get_data(datname);
|
||||
mesh_data_value_e value_type = data_ptr->get_valtype();
|
||||
mesh_data_type_e data_type = data_ptr->get_dattype();
|
||||
meshdata &data = get_data(datname);
|
||||
mesh_data_value_e value_type = data.valtype_;
|
||||
mesh_data_type_e data_type = data.loctype_;
|
||||
|
||||
if(value_type != Scalar)
|
||||
{
|
||||
@ -207,16 +204,13 @@ void gctl::gm_regular_grid::rtp(std::string datname, std::string rtpname, double
|
||||
}
|
||||
|
||||
// 检查是否存在与梯度数据同名的数据 若有则检查是否需要重新建立数据
|
||||
meshdata *new_data_ptr = add_data(rtpname, data_type, true, value_type);
|
||||
|
||||
array<double> *inval_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
array<double> *outval_ptr = (array<double>*) new_data_ptr->get_datval_ptr();
|
||||
meshdata &new_data = add_data(data_type, value_type, rtpname, 0.0);
|
||||
gctl::array<double> ingrid_ex;
|
||||
|
||||
int M = -1, N = -1;
|
||||
int ori_row, ori_col;
|
||||
if (data_type == NodeData) gctl::cosine_extrapolate_2d(*inval_ptr, ingrid_ex, rg_ynum, rg_xnum, M, N, ori_row, ori_col);
|
||||
else gctl::cosine_extrapolate_2d(*inval_ptr, ingrid_ex, rg_ynum-1, rg_xnum-1, M, N, ori_row, ori_col);
|
||||
if (data_type == NodeData) gctl::cosine_extrapolate_2d(data.datval_, ingrid_ex, rg_ynum, rg_xnum, M, N, ori_row, ori_col);
|
||||
else gctl::cosine_extrapolate_2d(data.datval_, ingrid_ex, rg_ynum-1, rg_xnum-1, M, N, ori_row, ori_col);
|
||||
|
||||
fftw_complex *in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*M*N);
|
||||
fftw_complex *out= (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*M*N);
|
||||
@ -312,7 +306,7 @@ void gctl::gm_regular_grid::rtp(std::string datname, std::string rtpname, double
|
||||
{
|
||||
for (int j = 0; j < rg_xnum; j++)
|
||||
{
|
||||
outval_ptr->at(i*rg_xnum+j) = ingrid_ex[(i+ori_row)*N+j+ori_col];
|
||||
new_data.datval_[i*rg_xnum+j] = ingrid_ex[(i+ori_row)*N+j+ori_col];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -322,7 +316,7 @@ void gctl::gm_regular_grid::rtp(std::string datname, std::string rtpname, double
|
||||
{
|
||||
for (int j = 0; j < rg_xnum-1; j++)
|
||||
{
|
||||
outval_ptr->at(i*(rg_xnum-1)+j) = ingrid_ex[(i+ori_row)*N+j+ori_col];
|
||||
new_data.datval_[i*(rg_xnum-1)+j] = ingrid_ex[(i+ori_row)*N+j+ori_col];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -337,17 +331,17 @@ void gctl::gm_regular_grid::drtp(std::string datname, std::string incname,
|
||||
std::string decname, std::string drtpname, int order)
|
||||
{
|
||||
// 首先获取数据并确定所有数据类型和格式都是一致的
|
||||
meshdata *md_ptr = get_data(datname);
|
||||
mesh_data_value_e dat_valtype = md_ptr->get_valtype();
|
||||
mesh_data_type_e dat_datype = md_ptr->get_dattype();
|
||||
meshdata *md_ptr = get_data_ptr(datname);
|
||||
mesh_data_value_e dat_valtype = md_ptr->valtype_;
|
||||
mesh_data_type_e dat_datype = md_ptr->loctype_;
|
||||
|
||||
md_ptr = get_data(incname);
|
||||
mesh_data_value_e inc_valtype = md_ptr->get_valtype();
|
||||
mesh_data_type_e inc_datype = md_ptr->get_dattype();
|
||||
md_ptr = get_data_ptr(incname);
|
||||
mesh_data_value_e inc_valtype = md_ptr->valtype_;
|
||||
mesh_data_type_e inc_datype = md_ptr->loctype_;
|
||||
|
||||
md_ptr = get_data(decname);
|
||||
mesh_data_value_e dec_valtype = md_ptr->get_valtype();
|
||||
mesh_data_type_e dec_datype = md_ptr->get_dattype();
|
||||
md_ptr = get_data_ptr(decname);
|
||||
mesh_data_value_e dec_valtype = md_ptr->valtype_;
|
||||
mesh_data_type_e dec_datype = md_ptr->loctype_;
|
||||
|
||||
std::string err_str;
|
||||
if (dat_valtype != Scalar)
|
||||
@ -363,32 +357,35 @@ void gctl::gm_regular_grid::drtp(std::string datname, std::string incname,
|
||||
}
|
||||
|
||||
// 获取数据指针
|
||||
md_ptr = get_data(datname);
|
||||
array<double> *dat_ptr = (array<double>*) md_ptr->get_datval_ptr();
|
||||
md_ptr = get_data_ptr(datname);
|
||||
array<double> *dat_ptr = &md_ptr->datval_;
|
||||
|
||||
md_ptr = get_data(incname);
|
||||
array<double> *inc_ptr = (array<double>*) md_ptr->get_datval_ptr();
|
||||
md_ptr = get_data_ptr(incname);
|
||||
array<double> *inc_ptr = &md_ptr->datval_;
|
||||
|
||||
md_ptr = get_data(decname);
|
||||
array<double> *dec_ptr = (array<double>*) md_ptr->get_datval_ptr();
|
||||
md_ptr = get_data_ptr(decname);
|
||||
array<double> *dec_ptr = &md_ptr->datval_;
|
||||
|
||||
// 检查是否存在同名的数据 若有则检查是否需要重新建立数据
|
||||
md_ptr = add_data(drtpname, dat_datype, true, dat_valtype);
|
||||
add_data(dat_datype, dat_valtype, drtpname, 0.0);
|
||||
md_ptr = get_data_ptr(drtpname);
|
||||
// 获取数据指针
|
||||
array<double> *drtp_ptr = (array<double>*) md_ptr->get_datval_ptr();
|
||||
array<double> *drtp_ptr = &md_ptr->datval_;
|
||||
|
||||
// 新建临时数据并获取数据指针
|
||||
std::string tmp_name = drtpname+"tmp";
|
||||
md_ptr = add_data(tmp_name, dat_datype, true, dat_valtype);
|
||||
array<double> *tmp_ptr = (array<double>*) md_ptr->get_datval_ptr();
|
||||
add_data(dat_datype, dat_valtype, tmp_name, 0.0);
|
||||
md_ptr = get_data_ptr(tmp_name);
|
||||
array<double> *tmp_ptr = &md_ptr->datval_;
|
||||
|
||||
std::string diff1_name = drtpname+"diff1";
|
||||
md_ptr = add_data(diff1_name, dat_datype, true, dat_valtype);
|
||||
array<double> *diff1_ptr = (array<double>*) md_ptr->get_datval_ptr();
|
||||
add_data(dat_datype, dat_valtype, diff1_name, 0.0);
|
||||
md_ptr = get_data_ptr(diff1_name);
|
||||
array<double> *diff1_ptr = &md_ptr->datval_;
|
||||
|
||||
std::string diff2_name = drtpname+"diff2";
|
||||
md_ptr = add_data(diff2_name, dat_datype, true, dat_valtype);
|
||||
array<double> *diff2_ptr = (array<double>*) md_ptr->get_datval_ptr();
|
||||
add_data(dat_datype, dat_valtype, diff2_name, 0.0);
|
||||
array<double> *diff2_ptr = &md_ptr->datval_;
|
||||
|
||||
// 计算平均磁化参数
|
||||
double meaninc = 0.0, meandec = 0.0;
|
||||
@ -461,9 +458,9 @@ void gctl::gm_regular_grid::drtp(std::string datname, std::string incname,
|
||||
|
||||
void gctl::gm_regular_grid::conti(std::string datname, std::string retname, double height)
|
||||
{
|
||||
meshdata *data_ptr = get_data(datname);
|
||||
mesh_data_value_e value_type = data_ptr->get_valtype();
|
||||
mesh_data_type_e data_type = data_ptr->get_dattype();
|
||||
meshdata &data = get_data(datname);
|
||||
mesh_data_value_e value_type = data.valtype_;
|
||||
mesh_data_type_e data_type = data.loctype_;
|
||||
|
||||
if(value_type != Scalar)
|
||||
{
|
||||
@ -471,16 +468,13 @@ void gctl::gm_regular_grid::conti(std::string datname, std::string retname, doub
|
||||
}
|
||||
|
||||
// 检查是否存在与梯度数据同名的数据 若有则检查是否需要重新建立数据
|
||||
meshdata *new_data_ptr = add_data(retname, data_type, true, value_type);
|
||||
|
||||
array<double> *inval_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
array<double> *outval_ptr = (array<double>*) new_data_ptr->get_datval_ptr();
|
||||
meshdata &new_data = add_data(data_type, value_type, retname, 0.0);
|
||||
gctl::array<double> ingrid_ex;
|
||||
|
||||
int M = -1, N = -1;
|
||||
int ori_row, ori_col;
|
||||
if (data_type == NodeData) gctl::cosine_extrapolate_2d(*inval_ptr, ingrid_ex, rg_ynum, rg_xnum, M, N, ori_row, ori_col);
|
||||
else gctl::cosine_extrapolate_2d(*inval_ptr, ingrid_ex, rg_ynum-1, rg_xnum-1, M, N, ori_row, ori_col);
|
||||
if (data_type == NodeData) gctl::cosine_extrapolate_2d(data.datval_, ingrid_ex, rg_ynum, rg_xnum, M, N, ori_row, ori_col);
|
||||
else gctl::cosine_extrapolate_2d(data.datval_, ingrid_ex, rg_ynum-1, rg_xnum-1, M, N, ori_row, ori_col);
|
||||
|
||||
fftw_complex *in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*M*N);
|
||||
fftw_complex *out= (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*M*N);
|
||||
@ -544,7 +538,7 @@ void gctl::gm_regular_grid::conti(std::string datname, std::string retname, doub
|
||||
{
|
||||
for (int j = 0; j < rg_xnum; j++)
|
||||
{
|
||||
outval_ptr->at(i*rg_xnum+j) = ingrid_ex[(i+ori_row)*N+j+ori_col];
|
||||
new_data.datval_[i*rg_xnum+j] = ingrid_ex[(i+ori_row)*N+j+ori_col];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -554,7 +548,7 @@ void gctl::gm_regular_grid::conti(std::string datname, std::string retname, doub
|
||||
{
|
||||
for (int j = 0; j < rg_xnum-1; j++)
|
||||
{
|
||||
outval_ptr->at(i*(rg_xnum-1)+j) = ingrid_ex[(i+ori_row)*N+j+ori_col];
|
||||
new_data.datval_[i*(rg_xnum-1)+j] = ingrid_ex[(i+ori_row)*N+j+ori_col];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -569,9 +563,9 @@ void gctl::gm_regular_grid::conti(std::string datname, std::string retname, doub
|
||||
|
||||
void gctl::gm_regular_grid::trend(std::string datname, std::string regname, std::string resname, int wx_size, int wy_size, int x_order, int y_order)
|
||||
{
|
||||
meshdata *data_ptr = get_data(datname);
|
||||
mesh_data_value_e value_type = data_ptr->get_valtype();
|
||||
mesh_data_type_e data_type = data_ptr->get_dattype();
|
||||
meshdata &data = get_data(datname);
|
||||
mesh_data_value_e value_type = data.valtype_;
|
||||
mesh_data_type_e data_type = data.loctype_;
|
||||
|
||||
if(value_type != Scalar)
|
||||
{
|
||||
@ -579,11 +573,8 @@ void gctl::gm_regular_grid::trend(std::string datname, std::string regname, std:
|
||||
}
|
||||
|
||||
// 检查是否存在与趋势场局部场数据同名的数据 若有则检查是否需要重新建立数据
|
||||
meshdata *reg_data_ptr = add_data(regname, data_type, true, value_type);
|
||||
meshdata *res_data_ptr = add_data(resname, data_type, true, value_type);
|
||||
array<double> *inval_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
array<double> *regval_ptr = (array<double>*) reg_data_ptr->get_datval_ptr();
|
||||
array<double> *resval_ptr = (array<double>*) res_data_ptr->get_datval_ptr();
|
||||
meshdata ®_data = add_data(data_type, value_type, regname, 0.0);
|
||||
meshdata &res_data = add_data(data_type, value_type, resname, 0.0);
|
||||
|
||||
matrix<double> data_mat;
|
||||
if (data_type == NodeData)
|
||||
@ -593,7 +584,7 @@ void gctl::gm_regular_grid::trend(std::string datname, std::string regname, std:
|
||||
{
|
||||
for (size_t j = 0; j < rg_xnum; j++)
|
||||
{
|
||||
data_mat[i][j] = inval_ptr->at(i*rg_xnum+j);
|
||||
data_mat[i][j] = data.datval_[i*rg_xnum+j];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -604,7 +595,7 @@ void gctl::gm_regular_grid::trend(std::string datname, std::string regname, std:
|
||||
{
|
||||
for (size_t j = 0; j < rg_xnum-1; j++)
|
||||
{
|
||||
data_mat[i][j] = inval_ptr->at(i*(rg_xnum-1)+j);
|
||||
data_mat[i][j] = data.datval_[i*(rg_xnum-1)+j];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -617,8 +608,8 @@ void gctl::gm_regular_grid::trend(std::string datname, std::string regname, std:
|
||||
{
|
||||
for (size_t j = 0; j < rg_xnum; j++)
|
||||
{
|
||||
regval_ptr->at(i*rg_xnum+j) = data_mat[i][j];
|
||||
resval_ptr->at(i*rg_xnum+j) = inval_ptr->at(i*rg_xnum+j) - data_mat[i][j];
|
||||
reg_data.datval_[i*rg_xnum+j] = data_mat[i][j];
|
||||
res_data.datval_[i*rg_xnum+j] = data.datval_[i*rg_xnum+j] - data_mat[i][j];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -628,8 +619,8 @@ void gctl::gm_regular_grid::trend(std::string datname, std::string regname, std:
|
||||
{
|
||||
for (size_t j = 0; j < rg_xnum-1; j++)
|
||||
{
|
||||
regval_ptr->at(i*(rg_xnum-1)+j) = data_mat[i][j];
|
||||
resval_ptr->at(i*(rg_xnum-1)+j) = inval_ptr->at(i*(rg_xnum-1)+j) - data_mat[i][j];
|
||||
reg_data.datval_[i*(rg_xnum-1)+j] = data_mat[i][j];
|
||||
res_data.datval_[i*(rg_xnum-1)+j] = data.datval_[i*(rg_xnum-1)+j] - data_mat[i][j];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -48,16 +48,15 @@ void gctl::grav_regular_mesh_2d::gobser(array<double> &out_data, std::string dat
|
||||
{
|
||||
array<double> rho(get_elenum(), 0.0);
|
||||
|
||||
meshdata *data_ptr = get_data(data_name);
|
||||
mesh_data_value_e vtype = data_ptr->get_valtype();
|
||||
meshdata &data = get_data(data_name);
|
||||
mesh_data_value_e vtype = data.valtype_;
|
||||
if (vtype == Scalar)
|
||||
{
|
||||
array<double> *val_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
for (int i = 0; i < val_ptr->size(); i++)
|
||||
for (int i = 0; i < data.datval_.size(); i++)
|
||||
{
|
||||
if (!std::isnan(val_ptr->at(i)))
|
||||
if (!std::isnan(data.datval_[i]))
|
||||
{
|
||||
rho[i] = val_ptr->at(i);
|
||||
rho[i] = data.datval_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -48,16 +48,15 @@ void gctl::gm_regular_mesh_3d::gobser(array<double> &out_data, std::string data_
|
||||
{
|
||||
array<double> rho(get_elenum(), 0.0);
|
||||
|
||||
meshdata *data_ptr = get_data(data_name);
|
||||
mesh_data_value_e vtype = data_ptr->get_valtype();
|
||||
meshdata &data = get_data(data_name);
|
||||
mesh_data_value_e vtype = data.valtype_;
|
||||
if (vtype == Scalar)
|
||||
{
|
||||
array<double> *val_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
for (int i = 0; i < val_ptr->size(); i++)
|
||||
for (int i = 0; i < data.datval_.size(); i++)
|
||||
{
|
||||
if (!std::isnan(val_ptr->at(i)))
|
||||
if (!std::isnan(data.datval_[i]))
|
||||
{
|
||||
rho[i] = val_ptr->at(i);
|
||||
rho[i] = data.datval_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -116,16 +115,15 @@ void gctl::gm_regular_mesh_3d::magobser(array<double> &out_data, std::string dat
|
||||
{
|
||||
array<double> sus(get_elenum(), 0.0);
|
||||
|
||||
meshdata *data_ptr = get_data(data_name);
|
||||
mesh_data_value_e vtype = data_ptr->get_valtype();
|
||||
meshdata &data = get_data(data_name);
|
||||
mesh_data_value_e vtype = data.valtype_;
|
||||
if (vtype == Scalar)
|
||||
{
|
||||
array<double> *val_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
for (int i = 0; i < val_ptr->size(); i++)
|
||||
for (int i = 0; i < data.datval_.size(); i++)
|
||||
{
|
||||
if (!std::isnan(val_ptr->at(i)))
|
||||
if (!std::isnan(data.datval_[i]))
|
||||
{
|
||||
sus[i] = val_ptr->at(i);
|
||||
sus[i] = data.datval_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -156,16 +154,15 @@ void gctl::gm_regular_mesh_3d::magobser(array<double> &out_data, std::string dat
|
||||
{
|
||||
array<double> sus(get_elenum(), 0.0);
|
||||
|
||||
meshdata *data_ptr = get_data(data_name);
|
||||
mesh_data_value_e vtype = data_ptr->get_valtype();
|
||||
meshdata &data = get_data(data_name);
|
||||
mesh_data_value_e vtype = data.valtype_;
|
||||
if (vtype == Scalar)
|
||||
{
|
||||
array<double> *val_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
for (int i = 0; i < val_ptr->size(); i++)
|
||||
for (int i = 0; i < data.datval_.size(); i++)
|
||||
{
|
||||
if (!std::isnan(val_ptr->at(i)))
|
||||
if (!std::isnan(data.datval_[i]))
|
||||
{
|
||||
sus[i] = val_ptr->at(i);
|
||||
sus[i] = data.datval_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -50,16 +50,15 @@ void gctl::gm_regular_mesh_sph_3d::gobser(array<double> &out_data, std::string d
|
||||
{
|
||||
array<double> rho(get_elenum(), 0.0);
|
||||
|
||||
meshdata *data_ptr = get_data(data_name);
|
||||
mesh_data_value_e vtype = data_ptr->get_valtype();
|
||||
meshdata &data = get_data(data_name);
|
||||
mesh_data_value_e vtype = data.valtype_;
|
||||
if (vtype == Scalar)
|
||||
{
|
||||
array<double> *val_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
for (int i = 0; i < val_ptr->size(); i++)
|
||||
for (int i = 0; i < data.datval_.size(); i++)
|
||||
{
|
||||
if (!std::isnan(val_ptr->at(i)))
|
||||
if (!std::isnan(data.datval_[i]))
|
||||
{
|
||||
rho[i] = val_ptr->at(i);
|
||||
rho[i] = data.datval_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -104,16 +103,15 @@ void gctl::gm_regular_mesh_sph_3d::magobser(array<double> &out_data, std::string
|
||||
{
|
||||
array<double> sus(get_elenum(), 0.0);
|
||||
|
||||
meshdata *data_ptr = get_data(data_name);
|
||||
mesh_data_value_e vtype = data_ptr->get_valtype();
|
||||
meshdata &data = get_data(data_name);
|
||||
mesh_data_value_e vtype = data.valtype_;
|
||||
if (vtype == Scalar)
|
||||
{
|
||||
array<double> *val_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
for (int i = 0; i < val_ptr->size(); i++)
|
||||
for (int i = 0; i < data.datval_.size(); i++)
|
||||
{
|
||||
if (!std::isnan(val_ptr->at(i)))
|
||||
if (!std::isnan(data.datval_[i]))
|
||||
{
|
||||
sus[i] = val_ptr->at(i);
|
||||
sus[i] = data.datval_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -66,16 +66,15 @@ void gctl::gm_tetrahedron_mesh::gobser(array<double> &out_obs, std::string data_
|
||||
{
|
||||
array<double> rho(get_elenum(), 0.0);
|
||||
|
||||
meshdata *data_ptr = get_data(data_name);
|
||||
mesh_data_value_e vtype = data_ptr->get_valtype();
|
||||
meshdata &data = get_data(data_name);
|
||||
mesh_data_value_e vtype = data.valtype_;
|
||||
if (vtype == Scalar)
|
||||
{
|
||||
array<double> *val_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
for (int i = 0; i < val_ptr->size(); i++)
|
||||
for (int i = 0; i < data.datval_.size(); i++)
|
||||
{
|
||||
if (!std::isnan(val_ptr->at(i)))
|
||||
if (!std::isnan(data.datval_[i]))
|
||||
{
|
||||
rho[i] = val_ptr->at(i);
|
||||
rho[i] = data.datval_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -101,16 +100,15 @@ void gctl::gm_tetrahedron_mesh::gobser(array<point3dc> &out_obs, std::string dat
|
||||
{
|
||||
array<double> rho(get_elenum(), 0.0);
|
||||
|
||||
meshdata *data_ptr = get_data(data_name);
|
||||
mesh_data_value_e vtype = data_ptr->get_valtype();
|
||||
meshdata &data = get_data(data_name);
|
||||
mesh_data_value_e vtype = data.valtype_;
|
||||
if (vtype == Scalar)
|
||||
{
|
||||
array<double> *val_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
for (int i = 0; i < val_ptr->size(); i++)
|
||||
for (int i = 0; i < data.datval_.size(); i++)
|
||||
{
|
||||
if (!std::isnan(val_ptr->at(i)))
|
||||
if (!std::isnan(data.datval_[i]))
|
||||
{
|
||||
rho[i] = val_ptr->at(i);
|
||||
rho[i] = data.datval_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -136,16 +134,15 @@ void gctl::gm_tetrahedron_mesh::gobser(array<tensor> &out_obs, std::string data_
|
||||
{
|
||||
array<double> rho(get_elenum(), 0.0);
|
||||
|
||||
meshdata *data_ptr = get_data(data_name);
|
||||
mesh_data_value_e vtype = data_ptr->get_valtype();
|
||||
meshdata &data = get_data(data_name);
|
||||
mesh_data_value_e vtype = data.valtype_;
|
||||
if (vtype == Scalar)
|
||||
{
|
||||
array<double> *val_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
for (int i = 0; i < val_ptr->size(); i++)
|
||||
for (int i = 0; i < data.datval_.size(); i++)
|
||||
{
|
||||
if (!std::isnan(val_ptr->at(i)))
|
||||
if (!std::isnan(data.datval_[i]))
|
||||
{
|
||||
rho[i] = val_ptr->at(i);
|
||||
rho[i] = data.datval_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -171,16 +168,15 @@ void gctl::gm_tetrahedron_mesh::gobser(array<double> &out_obs, std::string data_
|
||||
{
|
||||
array<double> rho(get_elenum(), 0.0);
|
||||
|
||||
meshdata *data_ptr = get_data(data_name);
|
||||
mesh_data_value_e vtype = data_ptr->get_valtype();
|
||||
meshdata &data = get_data(data_name);
|
||||
mesh_data_value_e vtype = data.valtype_;
|
||||
if (vtype == Scalar)
|
||||
{
|
||||
array<double> *val_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
for (int i = 0; i < val_ptr->size(); i++)
|
||||
for (int i = 0; i < data.datval_.size(); i++)
|
||||
{
|
||||
if (!std::isnan(val_ptr->at(i)))
|
||||
if (!std::isnan(data.datval_[i]))
|
||||
{
|
||||
rho[i] = val_ptr->at(i);
|
||||
rho[i] = data.datval_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -206,16 +202,15 @@ void gctl::gm_tetrahedron_mesh::gobser(array<point3dc> &out_obs, std::string dat
|
||||
{
|
||||
array<double> rho(get_elenum(), 0.0);
|
||||
|
||||
meshdata *data_ptr = get_data(data_name);
|
||||
mesh_data_value_e vtype = data_ptr->get_valtype();
|
||||
meshdata &data = get_data(data_name);
|
||||
mesh_data_value_e vtype = data.valtype_;
|
||||
if (vtype == Scalar)
|
||||
{
|
||||
array<double> *val_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
for (int i = 0; i < val_ptr->size(); i++)
|
||||
for (int i = 0; i < data.datval_.size(); i++)
|
||||
{
|
||||
if (!std::isnan(val_ptr->at(i)))
|
||||
if (!std::isnan(data.datval_[i]))
|
||||
{
|
||||
rho[i] = val_ptr->at(i);
|
||||
rho[i] = data.datval_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -241,16 +236,15 @@ void gctl::gm_tetrahedron_mesh::gobser(array<tensor> &out_obs, std::string data_
|
||||
{
|
||||
array<double> rho(get_elenum(), 0.0);
|
||||
|
||||
meshdata *data_ptr = get_data(data_name);
|
||||
mesh_data_value_e vtype = data_ptr->get_valtype();
|
||||
meshdata &data = get_data(data_name);
|
||||
mesh_data_value_e vtype = data.valtype_;
|
||||
if (vtype == Scalar)
|
||||
{
|
||||
array<double> *val_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
for (int i = 0; i < val_ptr->size(); i++)
|
||||
for (int i = 0; i < data.datval_.size(); i++)
|
||||
{
|
||||
if (!std::isnan(val_ptr->at(i)))
|
||||
if (!std::isnan(data.datval_[i]))
|
||||
{
|
||||
rho[i] = val_ptr->at(i);
|
||||
rho[i] = data.datval_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -275,19 +269,18 @@ void gctl::gm_tetrahedron_mesh::gobser(array<tensor> &out_obs, std::string data_
|
||||
void gctl::gm_tetrahedron_mesh::magkernel(matrix<double> &out_kernel, std::string sus_name, const array<point3dc> &obsp,
|
||||
double inclina_deg, double declina_deg, verbose_type_e verbose)
|
||||
{
|
||||
meshdata *data_ptr = get_data(sus_name);
|
||||
mesh_data_value_e vtype = data_ptr->get_valtype();
|
||||
meshdata &data = get_data(sus_name);
|
||||
mesh_data_value_e vtype = data.valtype_;
|
||||
|
||||
array<double> sus(get_elenum(), 0.0);
|
||||
|
||||
if (vtype == Scalar)
|
||||
{
|
||||
array<double> *val_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
for (int i = 0; i < val_ptr->size(); i++)
|
||||
for (int i = 0; i < data.datval_.size(); i++)
|
||||
{
|
||||
if (!std::isnan(val_ptr->at(i)))
|
||||
if (!std::isnan(data.datval_[i]))
|
||||
{
|
||||
sus[i] = val_ptr->at(i);
|
||||
sus[i] = data.datval_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -312,19 +305,18 @@ void gctl::gm_tetrahedron_mesh::magkernel(matrix<double> &out_kernel, std::strin
|
||||
void gctl::gm_tetrahedron_mesh::magkernel(matrix<point3dc> &out_kernel, std::string sus_name, const array<point3dc> &obsp,
|
||||
double inclina_deg, double declina_deg, verbose_type_e verbose)
|
||||
{
|
||||
meshdata *data_ptr = get_data(sus_name);
|
||||
mesh_data_value_e vtype = data_ptr->get_valtype();
|
||||
meshdata &data = get_data(sus_name);
|
||||
mesh_data_value_e vtype = data.valtype_;
|
||||
|
||||
array<double> sus(get_elenum(), 0.0);
|
||||
|
||||
if (vtype == Scalar)
|
||||
{
|
||||
array<double> *val_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
for (int i = 0; i < val_ptr->size(); i++)
|
||||
for (int i = 0; i < data.datval_.size(); i++)
|
||||
{
|
||||
if (!std::isnan(val_ptr->at(i)))
|
||||
if (!std::isnan(data.datval_[i]))
|
||||
{
|
||||
sus[i] = val_ptr->at(i);
|
||||
sus[i] = data.datval_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -349,19 +341,18 @@ void gctl::gm_tetrahedron_mesh::magkernel(matrix<point3dc> &out_kernel, std::str
|
||||
void gctl::gm_tetrahedron_mesh::magkernel(matrix<tensor> &out_kernel, std::string sus_name, const array<point3dc> &obsp,
|
||||
double inclina_deg, double declina_deg, verbose_type_e verbose)
|
||||
{
|
||||
meshdata *data_ptr = get_data(sus_name);
|
||||
mesh_data_value_e vtype = data_ptr->get_valtype();
|
||||
meshdata &data = get_data(sus_name);
|
||||
mesh_data_value_e vtype = data.valtype_;
|
||||
|
||||
array<double> sus(get_elenum(), 0.0);
|
||||
|
||||
if (vtype == Scalar)
|
||||
{
|
||||
array<double> *val_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
for (int i = 0; i < val_ptr->size(); i++)
|
||||
for (int i = 0; i < data.datval_.size(); i++)
|
||||
{
|
||||
if (!std::isnan(val_ptr->at(i)))
|
||||
if (!std::isnan(data.datval_[i]))
|
||||
{
|
||||
sus[i] = val_ptr->at(i);
|
||||
sus[i] = data.datval_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -386,19 +377,18 @@ void gctl::gm_tetrahedron_mesh::magkernel(matrix<tensor> &out_kernel, std::strin
|
||||
void gctl::gm_tetrahedron_mesh::magobser(array<double> &out_obs, std::string sus_name, const array<point3dc> &obsp,
|
||||
double inclina_deg, double declina_deg, verbose_type_e verbose)
|
||||
{
|
||||
meshdata *data_ptr = get_data(sus_name);
|
||||
mesh_data_value_e vtype = data_ptr->get_valtype();
|
||||
meshdata &data = get_data(sus_name);
|
||||
mesh_data_value_e vtype = data.valtype_;
|
||||
|
||||
array<double> sus(get_elenum(), 0.0);
|
||||
|
||||
if (vtype == Scalar)
|
||||
{
|
||||
array<double> *val_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
for (int i = 0; i < val_ptr->size(); i++)
|
||||
for (int i = 0; i < data.datval_.size(); i++)
|
||||
{
|
||||
if (!std::isnan(val_ptr->at(i)))
|
||||
if (!std::isnan(data.datval_[i]))
|
||||
{
|
||||
sus[i] = val_ptr->at(i);
|
||||
sus[i] = data.datval_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -423,19 +413,18 @@ void gctl::gm_tetrahedron_mesh::magobser(array<double> &out_obs, std::string sus
|
||||
void gctl::gm_tetrahedron_mesh::magobser(array<point3dc> &out_obs, std::string sus_name, const array<point3dc> &obsp,
|
||||
double inclina_deg, double declina_deg, verbose_type_e verbose)
|
||||
{
|
||||
meshdata *data_ptr = get_data(sus_name);
|
||||
mesh_data_value_e vtype = data_ptr->get_valtype();
|
||||
meshdata &data = get_data(sus_name);
|
||||
mesh_data_value_e vtype = data.valtype_;
|
||||
|
||||
array<double> sus(get_elenum(), 0.0);
|
||||
|
||||
if (vtype == Scalar)
|
||||
{
|
||||
array<double> *val_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
for (int i = 0; i < val_ptr->size(); i++)
|
||||
for (int i = 0; i < data.datval_.size(); i++)
|
||||
{
|
||||
if (!std::isnan(val_ptr->at(i)))
|
||||
if (!std::isnan(data.datval_[i]))
|
||||
{
|
||||
sus[i] = val_ptr->at(i);
|
||||
sus[i] = data.datval_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -460,19 +449,18 @@ void gctl::gm_tetrahedron_mesh::magobser(array<point3dc> &out_obs, std::string s
|
||||
void gctl::gm_tetrahedron_mesh::magobser(array<tensor> &out_obs, std::string sus_name, const array<point3dc> &obsp,
|
||||
double inclina_deg, double declina_deg, verbose_type_e verbose)
|
||||
{
|
||||
meshdata *data_ptr = get_data(sus_name);
|
||||
mesh_data_value_e vtype = data_ptr->get_valtype();
|
||||
meshdata &data = get_data(sus_name);
|
||||
mesh_data_value_e vtype = data.valtype_;
|
||||
|
||||
array<double> sus(get_elenum(), 0.0);
|
||||
|
||||
if (vtype == Scalar)
|
||||
{
|
||||
array<double> *val_ptr = (array<double>*) data_ptr->get_datval_ptr();
|
||||
for (int i = 0; i < val_ptr->size(); i++)
|
||||
for (int i = 0; i < data.datval_.size(); i++)
|
||||
{
|
||||
if (!std::isnan(val_ptr->at(i)))
|
||||
if (!std::isnan(data.datval_[i]))
|
||||
{
|
||||
sus[i] = val_ptr->at(i);
|
||||
sus[i] = data.datval_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -496,20 +484,12 @@ void gctl::gm_tetrahedron_mesh::magobser(array<tensor> &out_obs, std::string sus
|
||||
|
||||
void gctl::gm_tetrahedron_mesh::magobser(array<double> &out_obs, std::string magz_name, const array<point3dc> &obsp, verbose_type_e verbose)
|
||||
{
|
||||
meshdata *data_ptr;
|
||||
mesh_data_value_e vtype;
|
||||
|
||||
array<point3dc> magz(get_elenum(), point3dc(0.0, 0.0, 0.0));
|
||||
|
||||
data_ptr = get_data(magz_name);
|
||||
vtype = data_ptr->get_valtype();
|
||||
meshdata &data = get_data(magz_name);
|
||||
mesh_data_value_e vtype = data.valtype_;
|
||||
if (vtype == Vector)
|
||||
{
|
||||
array<point3dc> *val_ptr = (array<point3dc>*) data_ptr->get_datval_ptr();
|
||||
for (int i = 0; i < val_ptr->size(); i++)
|
||||
{
|
||||
magz[i] = val_ptr->at(i);
|
||||
}
|
||||
magz = data.export_vector();
|
||||
}
|
||||
else
|
||||
{
|
||||
@ -531,20 +511,12 @@ void gctl::gm_tetrahedron_mesh::magobser(array<double> &out_obs, std::string mag
|
||||
|
||||
void gctl::gm_tetrahedron_mesh::magobser(array<point3dc> &out_obs, std::string magz_name, const array<point3dc> &obsp, verbose_type_e verbose)
|
||||
{
|
||||
meshdata *data_ptr;
|
||||
mesh_data_value_e vtype;
|
||||
|
||||
array<point3dc> magz(get_elenum(), point3dc(0.0, 0.0, 0.0));
|
||||
|
||||
data_ptr = get_data(magz_name);
|
||||
vtype = data_ptr->get_valtype();
|
||||
meshdata &data = get_data(magz_name);
|
||||
mesh_data_value_e vtype = data.valtype_;
|
||||
if (vtype == Vector)
|
||||
{
|
||||
array<point3dc> *val_ptr = (array<point3dc>*) data_ptr->get_datval_ptr();
|
||||
for (int i = 0; i < val_ptr->size(); i++)
|
||||
{
|
||||
magz[i] = val_ptr->at(i);
|
||||
}
|
||||
magz = data.export_vector();
|
||||
}
|
||||
else
|
||||
{
|
||||
@ -566,20 +538,12 @@ void gctl::gm_tetrahedron_mesh::magobser(array<point3dc> &out_obs, std::string m
|
||||
|
||||
void gctl::gm_tetrahedron_mesh::magobser(array<tensor> &out_obs, std::string magz_name, const array<point3dc> &obsp, verbose_type_e verbose)
|
||||
{
|
||||
meshdata *data_ptr;
|
||||
mesh_data_value_e vtype;
|
||||
|
||||
array<point3dc> magz(get_elenum(), point3dc(0.0, 0.0, 0.0));
|
||||
|
||||
data_ptr = get_data(magz_name);
|
||||
vtype = data_ptr->get_valtype();
|
||||
meshdata &data = get_data(magz_name);
|
||||
mesh_data_value_e vtype = data.valtype_;
|
||||
if (vtype == Vector)
|
||||
{
|
||||
array<point3dc> *val_ptr = (array<point3dc>*) data_ptr->get_datval_ptr();
|
||||
for (int i = 0; i < val_ptr->size(); i++)
|
||||
{
|
||||
magz[i] = val_ptr->at(i);
|
||||
}
|
||||
magz = data.export_vector();
|
||||
}
|
||||
else
|
||||
{
|
||||
|
||||
@ -40,7 +40,7 @@ namespace gctl
|
||||
protected:
|
||||
array<grav_tetrahedron> grav_ele_;
|
||||
array<gravtet_para> grav_para_;
|
||||
array<mag_tetrahedron_ren17> mag_ele_;
|
||||
array<magtet_ren17> mag_ele_;
|
||||
array<magtet_para_ren17> mag_para_;
|
||||
|
||||
public:
|
||||
|
||||
@ -30,7 +30,7 @@
|
||||
|
||||
#define GCTL_MAG_TETRA_TOL 1e-16
|
||||
|
||||
void gctl::callink_magnetic_para_direct(array<mag_tetrahedron_ren17> &in_tet, array<magtet_para_ren17> &out_para, const array<point3dc> &magz)
|
||||
void gctl::callink_magnetic_para_direct(array<magtet_ren17> &in_tet, array<magtet_para_ren17> &out_para, const array<point3dc> &magz)
|
||||
{
|
||||
point3dc v1, v2, v3, ne, nf;
|
||||
|
||||
@ -66,7 +66,7 @@ void gctl::callink_magnetic_para_direct(array<mag_tetrahedron_ren17> &in_tet, ar
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::callink_magnetic_para_earth(array<mag_tetrahedron_ren17> &in_tet, array<magtet_para_ren17> &out_para,
|
||||
void gctl::callink_magnetic_para_earth(array<magtet_ren17> &in_tet, array<magtet_para_ren17> &out_para,
|
||||
double inclina_deg, double declina_deg, double field_tense)
|
||||
{
|
||||
if (declina_deg < -180.0 || declina_deg > 180.0 ||
|
||||
@ -116,7 +116,7 @@ void gctl::callink_magnetic_para_earth(array<mag_tetrahedron_ren17> &in_tet, arr
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::callink_magnetic_para_earth_sph(array<mag_tetrahedron_ren17> &in_tet, array<magtet_para_ren17> &out_para,
|
||||
void gctl::callink_magnetic_para_earth_sph(array<magtet_ren17> &in_tet, array<magtet_para_ren17> &out_para,
|
||||
double inclina_deg, double declina_deg, array<point3dc> *mag_vec, double field_tense)
|
||||
{
|
||||
if (declina_deg < -180.0 || declina_deg > 180.0 ||
|
||||
@ -177,7 +177,7 @@ void gctl::callink_magnetic_para_earth_sph(array<mag_tetrahedron_ren17> &in_tet,
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::callink_magnetic_para_earth_sph(mag_tetrahedron_ren17 &in_tet, magtet_para_ren17 &out_para,
|
||||
void gctl::callink_magnetic_para_earth_sph(magtet_ren17 &in_tet, magtet_para_ren17 &out_para,
|
||||
double inclina_deg, double declina_deg, point3dc *mag_vec, double field_tense)
|
||||
{
|
||||
if (declina_deg < -180.0 || declina_deg > 180.0 || inclina_deg < -90 || inclina_deg > 90 || field_tense < 0.0)
|
||||
@ -232,12 +232,12 @@ void gctl::callink_magnetic_para_earth_sph(mag_tetrahedron_ren17 &in_tet, magtet
|
||||
}
|
||||
|
||||
// declare algorithm for individual element and observation point
|
||||
double magkernel_tetrahedron_potential_sig(const gctl::mag_tetrahedron_ren17 &ele_ptr, const gctl::point3dc &op_ptr);
|
||||
gctl::point3dc magkernel_tetrahedron_gradient_sig(const gctl::mag_tetrahedron_ren17 &ele_ptr, const gctl::point3dc &op_ptr);
|
||||
gctl::tensor magkernel_tetrahedron_tensor_sig(const gctl::mag_tetrahedron_ren17 &ele_ptr, const gctl::point3dc &op_ptr);
|
||||
double magkernel_tetrahedron_potential_sig(const gctl::magtet_ren17 &ele_ptr, const gctl::point3dc &op_ptr);
|
||||
gctl::point3dc magkernel_tetrahedron_gradient_sig(const gctl::magtet_ren17 &ele_ptr, const gctl::point3dc &op_ptr);
|
||||
gctl::tensor magkernel_tetrahedron_tensor_sig(const gctl::magtet_ren17 &ele_ptr, const gctl::point3dc &op_ptr);
|
||||
|
||||
// define functions
|
||||
void gctl::magkernel(matrix<double> &out_kernel, const array<mag_tetrahedron_ren17> &ele,
|
||||
void gctl::magkernel(matrix<double> &out_kernel, const array<magtet_ren17> &ele,
|
||||
const array<point3dc> &obsp, verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
@ -261,7 +261,7 @@ void gctl::magkernel(matrix<double> &out_kernel, const array<mag_tetrahedron_ren
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::magkernel(matrix<point3dc> &out_kernel, const array<mag_tetrahedron_ren17> &ele,
|
||||
void gctl::magkernel(matrix<point3dc> &out_kernel, const array<magtet_ren17> &ele,
|
||||
const array<point3dc> &obsp, verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
@ -287,7 +287,7 @@ void gctl::magkernel(matrix<point3dc> &out_kernel, const array<mag_tetrahedron_r
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::magkernel(matrix<tensor> &out_kernel, const array<mag_tetrahedron_ren17> &ele,
|
||||
void gctl::magkernel(matrix<tensor> &out_kernel, const array<magtet_ren17> &ele,
|
||||
const array<point3dc> &obsp, verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
@ -311,7 +311,7 @@ void gctl::magkernel(matrix<tensor> &out_kernel, const array<mag_tetrahedron_ren
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::magkernel(matrix<double> &out_kernel, const array<mag_tetrahedron_ren17> &ele, const array<point3dc> &obsp,
|
||||
void gctl::magkernel(matrix<double> &out_kernel, const array<magtet_ren17> &ele, const array<point3dc> &obsp,
|
||||
double geo_declina, double geo_inclina, verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
@ -342,7 +342,7 @@ void gctl::magkernel(matrix<double> &out_kernel, const array<mag_tetrahedron_ren
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::magkernel(matrix<point3dc> &out_kernel, const array<mag_tetrahedron_ren17> &ele,
|
||||
void gctl::magkernel(matrix<point3dc> &out_kernel, const array<magtet_ren17> &ele,
|
||||
const array<point3ds> &obsp, verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
@ -391,7 +391,7 @@ void gctl::magkernel(matrix<point3dc> &out_kernel, const array<mag_tetrahedron_r
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::magkernel(matrix<tensor> &out_kernel, const array<mag_tetrahedron_ren17> &ele,
|
||||
void gctl::magkernel(matrix<tensor> &out_kernel, const array<magtet_ren17> &ele,
|
||||
const array<point3ds> &obsp, verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
@ -446,7 +446,7 @@ void gctl::magkernel(matrix<tensor> &out_kernel, const array<mag_tetrahedron_ren
|
||||
return;
|
||||
}
|
||||
|
||||
gctl::point3dc gctl::magkernel_single(const mag_tetrahedron_ren17 &ele, const point3ds &obsp)
|
||||
gctl::point3dc gctl::magkernel_single(const magtet_ren17 &ele, const point3ds &obsp)
|
||||
{
|
||||
tensor R;
|
||||
R[0][0] = sin((0.5-obsp.lat/180.0)*M_PI)*cos((2.0+obsp.lon/180.0)*M_PI);
|
||||
@ -469,7 +469,7 @@ gctl::point3dc gctl::magkernel_single(const mag_tetrahedron_ren17 &ele, const po
|
||||
return out_k;
|
||||
}
|
||||
|
||||
void gctl::magobser(array<double> &out_obs, const array<mag_tetrahedron_ren17> &ele, const array<point3dc> &obsp, verbose_type_e verbose)
|
||||
void gctl::magobser(array<double> &out_obs, const array<magtet_ren17> &ele, const array<point3dc> &obsp, verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
@ -492,7 +492,7 @@ void gctl::magobser(array<double> &out_obs, const array<mag_tetrahedron_ren17> &
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::magobser(array<point3dc> &out_obs, const array<mag_tetrahedron_ren17> &ele, const array<point3dc> &obsp, verbose_type_e verbose)
|
||||
void gctl::magobser(array<point3dc> &out_obs, const array<magtet_ren17> &ele, const array<point3dc> &obsp, verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
@ -517,7 +517,7 @@ void gctl::magobser(array<point3dc> &out_obs, const array<mag_tetrahedron_ren17>
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::magobser(array<tensor> &out_obs, const array<mag_tetrahedron_ren17> &ele, const array<point3dc> &obsp, verbose_type_e verbose)
|
||||
void gctl::magobser(array<tensor> &out_obs, const array<magtet_ren17> &ele, const array<point3dc> &obsp, verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
@ -540,7 +540,7 @@ void gctl::magobser(array<tensor> &out_obs, const array<mag_tetrahedron_ren17> &
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::magobser(array<double> &out_obs, const array<mag_tetrahedron_ren17> &ele, const array<point3dc> &obsp,
|
||||
void gctl::magobser(array<double> &out_obs, const array<magtet_ren17> &ele, const array<point3dc> &obsp,
|
||||
const array<double> &sus, verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
@ -564,7 +564,7 @@ void gctl::magobser(array<double> &out_obs, const array<mag_tetrahedron_ren17> &
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::magobser(array<point3dc> &out_obs, const array<mag_tetrahedron_ren17> &ele, const array<point3dc> &obsp,
|
||||
void gctl::magobser(array<point3dc> &out_obs, const array<magtet_ren17> &ele, const array<point3dc> &obsp,
|
||||
const array<double> &sus, verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
@ -590,7 +590,7 @@ void gctl::magobser(array<point3dc> &out_obs, const array<mag_tetrahedron_ren17>
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::magobser(array<tensor> &out_obs, const array<mag_tetrahedron_ren17> &ele, const array<point3dc> &obsp,
|
||||
void gctl::magobser(array<tensor> &out_obs, const array<magtet_ren17> &ele, const array<point3dc> &obsp,
|
||||
const array<double> &sus, verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
@ -614,7 +614,7 @@ void gctl::magobser(array<tensor> &out_obs, const array<mag_tetrahedron_ren17> &
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::magobser(array<point3dc> &out_obs, const array<mag_tetrahedron_ren17> &ele, const array<point3ds> &obsp,
|
||||
void gctl::magobser(array<point3dc> &out_obs, const array<magtet_ren17> &ele, const array<point3ds> &obsp,
|
||||
const array<double> &sus, verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
@ -660,7 +660,7 @@ void gctl::magobser(array<point3dc> &out_obs, const array<mag_tetrahedron_ren17>
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::magobser(array<tensor> &out_obs, const array<mag_tetrahedron_ren17> &ele, const array<point3ds> &obsp,
|
||||
void gctl::magobser(array<tensor> &out_obs, const array<magtet_ren17> &ele, const array<point3ds> &obsp,
|
||||
const array<double> &sus, verbose_type_e verbose)
|
||||
{
|
||||
int i, j;
|
||||
@ -712,7 +712,7 @@ void gctl::magobser(array<tensor> &out_obs, const array<mag_tetrahedron_ren17> &
|
||||
return;
|
||||
}
|
||||
|
||||
double magkernel_tetrahedron_potential_sig(const gctl::mag_tetrahedron_ren17 &tet, const gctl::point3dc &site)
|
||||
double magkernel_tetrahedron_potential_sig(const gctl::magtet_ren17 &tet, const gctl::point3dc &site)
|
||||
{
|
||||
double Rij_minus, Rij_plus, Sij_plus, Sij_minus, Rij0, mij0, wi0;
|
||||
double part1, part2, k0, absw, beta;
|
||||
@ -767,7 +767,7 @@ double magkernel_tetrahedron_potential_sig(const gctl::mag_tetrahedron_ren17 &te
|
||||
return out_pot;
|
||||
}
|
||||
|
||||
gctl::point3dc magkernel_tetrahedron_gradient_sig(const gctl::mag_tetrahedron_ren17 &tet, const gctl::point3dc &site)
|
||||
gctl::point3dc magkernel_tetrahedron_gradient_sig(const gctl::magtet_ren17 &tet, const gctl::point3dc &site)
|
||||
{
|
||||
double Rij_minus, Rij_plus, Sij_plus, Sij_minus, Rij0, mij0, wi0;
|
||||
double beta, Aij, sig, absw;
|
||||
@ -835,7 +835,7 @@ gctl::point3dc magkernel_tetrahedron_gradient_sig(const gctl::mag_tetrahedron_re
|
||||
return out_grad;
|
||||
}
|
||||
|
||||
gctl::tensor magkernel_tetrahedron_tensor_sig(const gctl::mag_tetrahedron_ren17 &tet, const gctl::point3dc &site)
|
||||
gctl::tensor magkernel_tetrahedron_tensor_sig(const gctl::magtet_ren17 &tet, const gctl::point3dc &site)
|
||||
{
|
||||
double Rij_minus, Rij_plus, Sij_plus, Sij_minus, Rij0, mij0, wi0;
|
||||
double beta, sig, absw;
|
||||
|
||||
@ -40,7 +40,7 @@ namespace gctl
|
||||
point3dc te[12]; // 四面体边切线矢量
|
||||
};
|
||||
|
||||
typedef type_tetrahedron<magtet_para_ren17> mag_tetrahedron_ren17; ///< 带magtet_para_ren17属性的四面体结构体
|
||||
typedef type_tetrahedron<magtet_para_ren17> magtet_ren17; ///< 带magtet_para_ren17属性的四面体结构体
|
||||
|
||||
/**
|
||||
* @brief Calculate the magnetic parameters of given tetrahedral elements.
|
||||
@ -52,7 +52,7 @@ namespace gctl
|
||||
* @param out_para Output parameters
|
||||
* @param magz Magnetic magnetizations
|
||||
*/
|
||||
void callink_magnetic_para_direct(array<mag_tetrahedron_ren17> &in_tet, array<magtet_para_ren17> &out_para, const array<point3dc> &magz);
|
||||
void callink_magnetic_para_direct(array<magtet_ren17> &in_tet, array<magtet_para_ren17> &out_para, const array<point3dc> &magz);
|
||||
|
||||
/**
|
||||
* @brief Calculate the magnetic parameters of given tetrahedral elements.
|
||||
@ -66,7 +66,7 @@ namespace gctl
|
||||
* @param declina_deg declination angle
|
||||
* @param field_tense Tense of the magnetic field
|
||||
*/
|
||||
void callink_magnetic_para_earth(array<mag_tetrahedron_ren17> &in_tet, array<magtet_para_ren17> &out_para,
|
||||
void callink_magnetic_para_earth(array<magtet_ren17> &in_tet, array<magtet_para_ren17> &out_para,
|
||||
double inclina_deg, double declina_deg, double field_tense = GCTL_T0);
|
||||
|
||||
/**
|
||||
@ -82,7 +82,7 @@ namespace gctl
|
||||
* @param mag_vec Output magnetization vectors (This is useful for data visualization)
|
||||
* @param field_tense Tense of the Earth's magnetic field
|
||||
*/
|
||||
void callink_magnetic_para_earth_sph(array<mag_tetrahedron_ren17> &in_tet, array<magtet_para_ren17> &out_para,
|
||||
void callink_magnetic_para_earth_sph(array<magtet_ren17> &in_tet, array<magtet_para_ren17> &out_para,
|
||||
double inclina_deg, double declina_deg, array<point3dc> *mag_vec = nullptr, double field_tense = GCTL_T0);
|
||||
|
||||
/**
|
||||
@ -98,52 +98,52 @@ namespace gctl
|
||||
* @param mag_vec Output magnetization vectors (This is useful for data visualization)
|
||||
* @param field_tense Tense of the Earth's magnetic field
|
||||
*/
|
||||
void callink_magnetic_para_earth_sph(mag_tetrahedron_ren17 &in_tet, magtet_para_ren17 &out_para,
|
||||
void callink_magnetic_para_earth_sph(magtet_ren17 &in_tet, magtet_para_ren17 &out_para,
|
||||
double inclina_deg, double declina_deg, point3dc *mag_vec = nullptr, double field_tense = GCTL_T0);
|
||||
|
||||
/**
|
||||
* @brief Calculate the kernel matrix of the magnetic poential data.
|
||||
*
|
||||
* @note Use callink_magnetic_para_earth() to initialize the mag_tetrahedron_ren17 elements.
|
||||
* @note Use callink_magnetic_para_earth() to initialize the magtet_ren17 elements.
|
||||
*
|
||||
* @param out_kernel Output kernel matrix
|
||||
* @param ele Magnetic tetrahedrons.
|
||||
* @param obsp Observation sites
|
||||
* @param verbose Output info level
|
||||
*/
|
||||
void magkernel(matrix<double> &out_kernel, const array<mag_tetrahedron_ren17> &ele, const array<point3dc> &obsp,
|
||||
void magkernel(matrix<double> &out_kernel, const array<magtet_ren17> &ele, const array<point3dc> &obsp,
|
||||
verbose_type_e verbose = FullMsg);
|
||||
|
||||
/**
|
||||
* @brief Calculate the kernel matrix of the magnetic conponents data.
|
||||
*
|
||||
* @note Use callink_magnetic_para_earth() to initialize the mag_tetrahedron_ren17 elements.
|
||||
* @note Use callink_magnetic_para_earth() to initialize the magtet_ren17 elements.
|
||||
*
|
||||
* @param out_kernel Output kernel matrix. Directional components are stored accordingly.
|
||||
* @param ele Magnetic tetrahedrons.
|
||||
* @param obsp Observation sites
|
||||
* @param verbose Output info level
|
||||
*/
|
||||
void magkernel(matrix<point3dc> &out_kernel, const array<mag_tetrahedron_ren17> &ele, const array<point3dc> &obsp,
|
||||
void magkernel(matrix<point3dc> &out_kernel, const array<magtet_ren17> &ele, const array<point3dc> &obsp,
|
||||
verbose_type_e verbose = FullMsg);
|
||||
|
||||
/**
|
||||
* @brief Calculate the kernel matrix of the magnetic tensor data.
|
||||
*
|
||||
* @note Use callink_magnetic_para_earth() to initialize the mag_tetrahedron_ren17 elements.
|
||||
* @note Use callink_magnetic_para_earth() to initialize the magtet_ren17 elements.
|
||||
*
|
||||
* @param out_kernel Output kernel matrix. Directional components are stored accordingly.
|
||||
* @param ele Magnetic tetrahedrons.
|
||||
* @param obsp Observation sites
|
||||
* @param verbose Output info level
|
||||
*/
|
||||
void magkernel(matrix<tensor> &out_kernel, const array<mag_tetrahedron_ren17> &ele, const array<point3dc> &obsp,
|
||||
void magkernel(matrix<tensor> &out_kernel, const array<magtet_ren17> &ele, const array<point3dc> &obsp,
|
||||
verbose_type_e verbose = FullMsg);
|
||||
|
||||
/**
|
||||
* @brief Calculate the kernel matrix of the deltaT anomaly data.
|
||||
*
|
||||
* @note Use callink_magnetic_para_earth() to initialize the mag_tetrahedron_ren17 elements.
|
||||
* @note Use callink_magnetic_para_earth() to initialize the magtet_ren17 elements.
|
||||
*
|
||||
* @param out_kernel Output kernel matrix
|
||||
* @param ele Magnetic tetrahedrons.
|
||||
@ -152,33 +152,33 @@ namespace gctl
|
||||
* @param geo_inclina Inclination angle of the geo-magnetic field.
|
||||
* @param verbose Output info level
|
||||
*/
|
||||
void magkernel(matrix<double> &out_kernel, const array<mag_tetrahedron_ren17> &ele, const array<point3dc> &obsp,
|
||||
void magkernel(matrix<double> &out_kernel, const array<magtet_ren17> &ele, const array<point3dc> &obsp,
|
||||
double geo_declina, double geo_inclina, verbose_type_e verbose = FullMsg);
|
||||
|
||||
/**
|
||||
* @brief Calculate the kernel matrix of the magnetic conponents data.
|
||||
*
|
||||
* @note Use callink_magnetic_para_earth_sph() to initialize the mag_tetrahedron_ren17 elements.
|
||||
* @note Use callink_magnetic_para_earth_sph() to initialize the magtet_ren17 elements.
|
||||
*
|
||||
* @param out_kernel Output kernel matrix. Directional components are stored accordingly.
|
||||
* @param ele Magnetic tetrahedrons.
|
||||
* @param obsp Observation sites
|
||||
* @param verbose Output info level
|
||||
*/
|
||||
void magkernel(matrix<point3dc> &out_kernel, const array<mag_tetrahedron_ren17> &ele, const array<point3ds> &obsp,
|
||||
void magkernel(matrix<point3dc> &out_kernel, const array<magtet_ren17> &ele, const array<point3ds> &obsp,
|
||||
verbose_type_e verbose = FullMsg);
|
||||
|
||||
/**
|
||||
* @brief Calculate the kernel matrix of the magnetic tensor data.
|
||||
*
|
||||
* @note Use callink_magnetic_para_earth_sph() to initialize the mag_tetrahedron_ren17 elements.
|
||||
* @note Use callink_magnetic_para_earth_sph() to initialize the magtet_ren17 elements.
|
||||
*
|
||||
* @param out_kernel Output kernel matrix. Directional components are stored accordingly.
|
||||
* @param ele Magnetic tetrahedrons.
|
||||
* @param obsp Observation sites
|
||||
* @param verbose Output info level
|
||||
*/
|
||||
void magkernel(matrix<tensor> &out_kernel, const array<mag_tetrahedron_ren17> &ele, const array<point3ds> &obsp,
|
||||
void magkernel(matrix<tensor> &out_kernel, const array<magtet_ren17> &ele, const array<point3ds> &obsp,
|
||||
verbose_type_e verbose = FullMsg);
|
||||
|
||||
/**
|
||||
@ -188,48 +188,48 @@ namespace gctl
|
||||
* @param obsp The observation point
|
||||
* @return point3dc The returned magnetic componments
|
||||
*/
|
||||
point3dc magkernel_single(const mag_tetrahedron_ren17 &ele, const point3ds &obsp);
|
||||
point3dc magkernel_single(const magtet_ren17 &ele, const point3ds &obsp);
|
||||
|
||||
/**
|
||||
* @brief Calculate magnetic potential data.
|
||||
*
|
||||
* @note The magnetic susceptibility is impliclity setted using the magnetic magnetization values. Use callink_magnetic_para_direct() to initialize the mag_tetrahedron_ren17 elements.
|
||||
* @note The magnetic susceptibility is impliclity setted using the magnetic magnetization values. Use callink_magnetic_para_direct() to initialize the magtet_ren17 elements.
|
||||
*
|
||||
* @param out_obs Output magnetic field data
|
||||
* @param ele Magnetic tetrahedrons.
|
||||
* @param obsp Observation sites
|
||||
* @param verbose Output info level
|
||||
*/
|
||||
void magobser(array<double> &out_obs, const array<mag_tetrahedron_ren17> &ele, const array<point3dc> &obsp, verbose_type_e verbose = FullMsg);
|
||||
void magobser(array<double> &out_obs, const array<magtet_ren17> &ele, const array<point3dc> &obsp, verbose_type_e verbose = FullMsg);
|
||||
|
||||
/**
|
||||
* @brief Calculate magnetic component data.
|
||||
*
|
||||
* @note The magnetic susceptibility is impliclity setted using the magnetic magnetization values. Use callink_magnetic_para_direct() to initialize the mag_tetrahedron_ren17 elements.
|
||||
* @note The magnetic susceptibility is impliclity setted using the magnetic magnetization values. Use callink_magnetic_para_direct() to initialize the magtet_ren17 elements.
|
||||
*
|
||||
* @param out_obs Output magnetic field data. Directional components are stored accordingly. Delta_T anomalies could be retrived using the function magnetic_components2deltaT().
|
||||
* @param ele Magnetic tetrahedrons.
|
||||
* @param obsp Observation sites
|
||||
* @param verbose Output info level
|
||||
*/
|
||||
void magobser(array<point3dc> &out_obs, const array<mag_tetrahedron_ren17> &ele, const array<point3dc> &obsp, verbose_type_e verbose = FullMsg);
|
||||
void magobser(array<point3dc> &out_obs, const array<magtet_ren17> &ele, const array<point3dc> &obsp, verbose_type_e verbose = FullMsg);
|
||||
|
||||
/**
|
||||
* @brief Calculate magnetic tensor data.
|
||||
*
|
||||
* @note The magnetic susceptibility is impliclity setted using the magnetic magnetization values. Use callink_magnetic_para_direct() to initialize the mag_tetrahedron_ren17 elements.
|
||||
* @note The magnetic susceptibility is impliclity setted using the magnetic magnetization values. Use callink_magnetic_para_direct() to initialize the magtet_ren17 elements.
|
||||
*
|
||||
* @param out_obs Output magnetic field data. Directional components are stored accordingly. Delta_T gradient anomalies could be retrived using the function magnetic_tensors2deltaTs().
|
||||
* @param ele Magnetic tetrahedrons.
|
||||
* @param obsp Observation sites
|
||||
* @param verbose Output info level
|
||||
*/
|
||||
void magobser(array<tensor> &out_obs, const array<mag_tetrahedron_ren17> &ele, const array<point3dc> &obsp, verbose_type_e verbose = FullMsg);
|
||||
void magobser(array<tensor> &out_obs, const array<magtet_ren17> &ele, const array<point3dc> &obsp, verbose_type_e verbose = FullMsg);
|
||||
|
||||
/**
|
||||
* @brief Calculate magnetic potential data.
|
||||
*
|
||||
* @note Use callink_magnetic_para_earth() to initialize the mag_tetrahedron_ren17 elements.
|
||||
* @note Use callink_magnetic_para_earth() to initialize the magtet_ren17 elements.
|
||||
*
|
||||
* @param out_obs Output magnetic field data
|
||||
* @param ele Magnetic tetrahedrons.
|
||||
@ -237,13 +237,13 @@ namespace gctl
|
||||
* @param sus Magnetic susceptibilities
|
||||
* @param verbose Output info level
|
||||
*/
|
||||
void magobser(array<double> &out_obs, const array<mag_tetrahedron_ren17> &ele, const array<point3dc> &obsp,
|
||||
void magobser(array<double> &out_obs, const array<magtet_ren17> &ele, const array<point3dc> &obsp,
|
||||
const array<double> &sus, verbose_type_e verbose = FullMsg);
|
||||
|
||||
/**
|
||||
* @brief Calculate magnetic component data.
|
||||
*
|
||||
* @note Use callink_magnetic_para_earth() to initialize the mag_tetrahedron_ren17 elements.
|
||||
* @note Use callink_magnetic_para_earth() to initialize the magtet_ren17 elements.
|
||||
*
|
||||
* @param out_obs Output magnetic field data. Directional components are stored accordingly. Delta_T anomalies could be retrived using the function magnetic_components2deltaT().
|
||||
* @param ele Magnetic tetrahedrons.
|
||||
@ -251,13 +251,13 @@ namespace gctl
|
||||
* @param sus Magnetic susceptibilities
|
||||
* @param verbose Output info level
|
||||
*/
|
||||
void magobser(array<point3dc> &out_obs, const array<mag_tetrahedron_ren17> &ele, const array<point3dc> &obsp,
|
||||
void magobser(array<point3dc> &out_obs, const array<magtet_ren17> &ele, const array<point3dc> &obsp,
|
||||
const array<double> &sus, verbose_type_e verbose = FullMsg);
|
||||
|
||||
/**
|
||||
* @brief Calculate magnetic tensor data.
|
||||
*
|
||||
* @note Use callink_magnetic_para_earth() to initialize the mag_tetrahedron_ren17 elements.
|
||||
* @note Use callink_magnetic_para_earth() to initialize the magtet_ren17 elements.
|
||||
*
|
||||
* @param out_obs Output magnetic field data. Directional tensors are stored accordingly. Delta_T gradient anomalies could be retrived using the function magnetic_tensors2deltaTs().
|
||||
* @param ele Magnetic tetrahedrons.
|
||||
@ -265,13 +265,13 @@ namespace gctl
|
||||
* @param sus Magnetic susceptibilities
|
||||
* @param verbose Output info level
|
||||
*/
|
||||
void magobser(array<tensor> &out_obs, const array<mag_tetrahedron_ren17> &ele, const array<point3dc> &obsp,
|
||||
void magobser(array<tensor> &out_obs, const array<magtet_ren17> &ele, const array<point3dc> &obsp,
|
||||
const array<double> &sus, verbose_type_e verbose = FullMsg);
|
||||
|
||||
/**
|
||||
* @brief Calculate magnetic component data.
|
||||
*
|
||||
* @note Use callink_magnetic_para_earth_sph() to initialize the mag_tetrahedron_ren17 elements.
|
||||
* @note Use callink_magnetic_para_earth_sph() to initialize the magtet_ren17 elements.
|
||||
*
|
||||
* @param out_obs Output magnetic field data. Directional components are stored accordingly. Delta_T anomalies could be retrived using the function magnetic_components2deltaT().
|
||||
* @param ele Magnetic tetrahedrons.
|
||||
@ -279,13 +279,13 @@ namespace gctl
|
||||
* @param sus Magnetic susceptibilities
|
||||
* @param verbose Output info level
|
||||
*/
|
||||
void magobser(array<point3dc> &out_obs, const array<mag_tetrahedron_ren17> &ele, const array<point3ds> &obsp,
|
||||
void magobser(array<point3dc> &out_obs, const array<magtet_ren17> &ele, const array<point3ds> &obsp,
|
||||
const array<double> &sus, verbose_type_e verbose = FullMsg);
|
||||
|
||||
/**
|
||||
* @brief Calculate magnetic tensor data.
|
||||
*
|
||||
* @note Use callink_magnetic_para_earth_sph() to initialize the mag_tetrahedron_ren17 elements.
|
||||
* @note Use callink_magnetic_para_earth_sph() to initialize the magtet_ren17 elements.
|
||||
*
|
||||
* @param out_obs Output magnetic field data. Directional tensors are stored accordingly. Delta_T gradient anomalies could be retrived using the function magnetic_tensors2deltaTs().
|
||||
* @param ele Magnetic tetrahedrons.
|
||||
@ -293,7 +293,7 @@ namespace gctl
|
||||
* @param sus Magnetic susceptibilities
|
||||
* @param verbose Output info level
|
||||
*/
|
||||
void magobser(array<tensor> &out_obs, const array<mag_tetrahedron_ren17> &ele, const array<point3ds> &obsp,
|
||||
void magobser(array<tensor> &out_obs, const array<magtet_ren17> &ele, const array<point3ds> &obsp,
|
||||
const array<double> &sus, verbose_type_e verbose = FullMsg);
|
||||
}
|
||||
#endif // _GCTL_MAG_KERNEL_TETRAHEDRON_REN2017_H
|
||||
@ -27,7 +27,7 @@
|
||||
|
||||
#include "mkernel_tricone.h"
|
||||
|
||||
void gctl::callink_magnetic_para(array<mag_tricone> &in_cone, array<magcone_para> &out_para, const array<point3dc> &mag_B)
|
||||
void gctl::callink_magnetic_para(array<magcone> &in_cone, array<magcone_para> &out_para, const array<point3dc> &mag_B)
|
||||
{
|
||||
point3dc v1, v2, v3, ne, nf, mag_z;
|
||||
|
||||
@ -63,7 +63,7 @@ void gctl::callink_magnetic_para(array<mag_tricone> &in_cone, array<magcone_para
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::callink_magnetic_para_earth_sph(array<mag_tricone> &in_tet, array<magcone_para> &out_para,
|
||||
void gctl::callink_magnetic_para_earth_sph(array<magcone> &in_tet, array<magcone_para> &out_para,
|
||||
double inclina_deg, double declina_deg, array<point3dc> *mag_vec, double field_tense)
|
||||
{
|
||||
if (declina_deg < -180.0 || declina_deg > 180.0 ||
|
||||
@ -122,7 +122,7 @@ void gctl::callink_magnetic_para_earth_sph(array<mag_tricone> &in_tet, array<mag
|
||||
return;
|
||||
}
|
||||
|
||||
gctl::point3dc gctl::magkernel_single(const mag_tricone &a_ele, const point3ds &a_op, tensor *R_ptr)
|
||||
gctl::point3dc gctl::magkernel_single(const magcone &a_ele, const point3ds &a_op, tensor *R_ptr)
|
||||
{
|
||||
int f,e;
|
||||
double Le,wf;
|
||||
@ -175,8 +175,8 @@ gctl::point3dc gctl::magkernel_single(const mag_tricone &a_ele, const point3ds &
|
||||
return out_c;
|
||||
}
|
||||
|
||||
void gctl::magkernel(matrix<double> &kernel, const array<mag_tricone> &top_ele,
|
||||
const array<mag_tricone> &btm_ele, const array<point3ds> &obsp,
|
||||
void gctl::magkernel(matrix<double> &kernel, const array<magcone> &top_ele,
|
||||
const array<magcone> &btm_ele, const array<point3ds> &obsp,
|
||||
magnetic_field_type_e comp_type, verbose_type_e verbose)
|
||||
{
|
||||
if (comp_type != Bx && comp_type != By && comp_type != Bz)
|
||||
@ -237,7 +237,7 @@ void gctl::magkernel(matrix<double> &kernel, const array<mag_tricone> &top_ele,
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::magkernel(spmat<double> &kernel, const array<mag_tricone> &top_ele, const array<mag_tricone> &btm_ele,
|
||||
void gctl::magkernel(spmat<double> &kernel, const array<magcone> &top_ele, const array<magcone> &btm_ele,
|
||||
const array<point3ds> &obsp, double cut_angle, magnetic_field_type_e comp_type, verbose_type_e verbose)
|
||||
{
|
||||
if (comp_type != Bx && comp_type != By && comp_type != Bz)
|
||||
@ -333,8 +333,8 @@ void gctl::magkernel(spmat<double> &kernel, const array<mag_tricone> &top_ele, c
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::magobser(array<point3dc> &out_obs, const array<mag_tricone> &top_ele,
|
||||
const array<mag_tricone> &btm_ele, const array<point3ds> &obsp, const array<double> &sus,
|
||||
void gctl::magobser(array<point3dc> &out_obs, const array<magcone> &top_ele,
|
||||
const array<magcone> &btm_ele, const array<point3ds> &obsp, const array<double> &sus,
|
||||
verbose_type_e verbose)
|
||||
{
|
||||
if (top_ele.size() != btm_ele.size())
|
||||
|
||||
@ -40,7 +40,7 @@ namespace gctl
|
||||
double edglen[12]; ///< edge lengths of six edges
|
||||
};
|
||||
|
||||
typedef type_tricone<magcone_para> mag_tricone; ///< 带magcone_para属性的三角锥结构体
|
||||
typedef type_tricone<magcone_para> magcone; ///< 带magcone_para属性的三角锥结构体
|
||||
|
||||
/**
|
||||
* @brief Calculate the magnetic parameters of given tricone elements.
|
||||
@ -49,7 +49,7 @@ namespace gctl
|
||||
* @param out_para Output parameters
|
||||
* @param mag_B magnetization vecrtors
|
||||
*/
|
||||
void callink_magnetic_para(array<mag_tricone> &in_cone, array<magcone_para> &out_para, const array<point3dc> &mag_B);
|
||||
void callink_magnetic_para(array<magcone> &in_cone, array<magcone_para> &out_para, const array<point3dc> &mag_B);
|
||||
|
||||
/**
|
||||
* @brief Calculate the magnetic parameters of given tricone elements wrt. the spherical coordinates.
|
||||
@ -64,18 +64,18 @@ namespace gctl
|
||||
* @param mag_vec Output magnetization vectors (This is useful for data visualization)
|
||||
* @param field_tense Tense of the Earth's magnetic field
|
||||
*/
|
||||
void callink_magnetic_para_earth_sph(array<mag_tricone> &in_tet, array<magcone_para> &out_para,
|
||||
void callink_magnetic_para_earth_sph(array<magcone> &in_tet, array<magcone_para> &out_para,
|
||||
double inclina_deg, double declina_deg, array<point3dc> *mag_vec = nullptr, double field_tense = GCTL_T0);
|
||||
|
||||
point3dc magkernel_single(const mag_tricone &a_ele, const point3ds &a_op, tensor *R_ptr = nullptr);
|
||||
point3dc magkernel_single(const magcone &a_ele, const point3ds &a_op, tensor *R_ptr = nullptr);
|
||||
|
||||
void magkernel(matrix<double> &kernel, const array<mag_tricone> &top_ele, const array<mag_tricone> &btm_ele,
|
||||
void magkernel(matrix<double> &kernel, const array<magcone> &top_ele, const array<magcone> &btm_ele,
|
||||
const array<point3ds> &obsp, magnetic_field_type_e comp_type = Bz, verbose_type_e verbose = FullMsg);
|
||||
|
||||
void magkernel(spmat<double> &kernel, const array<mag_tricone> &top_ele, const array<mag_tricone> &btm_ele,
|
||||
void magkernel(spmat<double> &kernel, const array<magcone> &top_ele, const array<magcone> &btm_ele,
|
||||
const array<point3ds> &obsp, double cut_angle, magnetic_field_type_e comp_type = Bz, verbose_type_e verbose = FullMsg);
|
||||
|
||||
void magobser(array<point3dc> &out_obs, const array<mag_tricone> &top_ele, const array<mag_tricone> &btm_ele,
|
||||
void magobser(array<point3dc> &out_obs, const array<magcone> &top_ele, const array<magcone> &btm_ele,
|
||||
const array<point3ds> &obsp, const array<double> &sus, verbose_type_e verbose = FullMsg);
|
||||
}
|
||||
|
||||
|
||||
589
lib/potential/mkernel_tricone_Ren2017.cpp
Normal file
589
lib/potential/mkernel_tricone_Ren2017.cpp
Normal file
@ -0,0 +1,589 @@
|
||||
/********************************************************
|
||||
* ██████╗ ██████╗████████╗██╗
|
||||
* ██╔════╝ ██╔════╝╚══██╔══╝██║
|
||||
* ██║ ███╗██║ ██║ ██║
|
||||
* ██║ ██║██║ ██║ ██║
|
||||
* ╚██████╔╝╚██████╗ ██║ ███████╗
|
||||
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
|
||||
* Geophysical Computational Tools & Library (GCTL)
|
||||
*
|
||||
* Copyright (c) 2022 Yi Zhang (yizhang-geo@zju.edu.cn)
|
||||
*
|
||||
* GCTL is distributed under a dual licensing scheme. You can redistribute
|
||||
* it and/or modify it under the terms of the GNU Lesser General Public
|
||||
* License as published by the Free Software Foundation, either version 2
|
||||
* of the License, or (at your option) any later version. You should have
|
||||
* received a copy of the GNU Lesser General Public License along with this
|
||||
* program. If not, see <http://www.gnu.org/licenses/>.
|
||||
*
|
||||
* If the terms and conditions of the LGPL v.2. would prevent you from using
|
||||
* the GCTL, please consider the option to obtain a commercial license for a
|
||||
* fee. These licenses are offered by the GCTL's original author. As a rule,
|
||||
* licenses are provided "as-is", unlimited in time for a one time fee. Please
|
||||
* send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget
|
||||
* to include some description of your company and the realm of its activities.
|
||||
* Also add information on how to contact you by electronic and paper mail.
|
||||
******************************************************/
|
||||
|
||||
#include "mkernel_tricone_Ren2017.h"
|
||||
#include "cmath"
|
||||
|
||||
double mag_tolerance = 1e-16;
|
||||
|
||||
void gctl::set_magcone_ren17_tolerance(double tol)
|
||||
{
|
||||
if (tol > 0) mag_tolerance = tol;
|
||||
else throw invalid_argument("Invalid tolerance value. From gctl::set_magcone_ren17_tolerance(...)");
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::callink_magnetic_para(array<magcone_ren17> &in_cone, array<magcone_para_ren17> &out_para, const array<point3dc> &mag_B)
|
||||
{
|
||||
point3dc v1, v2, v3, ne, nf;
|
||||
|
||||
out_para.resize(in_cone.size());
|
||||
for (int i = 0; i < in_cone.size(); ++i)
|
||||
{
|
||||
if (in_cone[i].vert[0] == nullptr || in_cone[i].vert[1] == nullptr ||
|
||||
in_cone[i].vert[2] == nullptr || in_cone[i].vert[3] == nullptr)
|
||||
{
|
||||
throw runtime_error("Invalid vertex pointer. From callink_magnetic_para(...)");
|
||||
}
|
||||
|
||||
for (int f = 0; f < 4; ++f)
|
||||
{
|
||||
v1 = *in_cone[i].fget(f, 1) - *in_cone[i].fget(f, 0);
|
||||
v2 = *in_cone[i].fget(f, 2) - *in_cone[i].fget(f, 0);
|
||||
nf = cross(v1, v2).normal();
|
||||
out_para[i].nf[f] = nf;
|
||||
|
||||
for (int e = 0; e < 3; ++e)
|
||||
{
|
||||
v3 = *in_cone[i].fget(f, (e+1)%3) - *in_cone[i].fget(f, e);
|
||||
ne = cross(v3, nf).normal();
|
||||
out_para[i].ne[e+f*3] = ne;
|
||||
out_para[i].te[e+f*3] = cross(nf, ne).normal();
|
||||
}
|
||||
|
||||
out_para[i].mag_amp[f] = dot(mag_B[i], nf);
|
||||
}
|
||||
|
||||
in_cone[i].att = out_para.get(i);
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::callink_magnetic_para_earth_sph(array<magcone_ren17> &in_cone, array<magcone_para_ren17> &out_para,
|
||||
double inclina_deg, double declina_deg, array<point3dc> *mag_vec, double field_tense)
|
||||
{
|
||||
if (declina_deg < -180.0 || declina_deg > 180.0 ||
|
||||
inclina_deg < -90 || inclina_deg > 90 || field_tense < 0.0)
|
||||
{
|
||||
throw invalid_argument("Invalid parameters. From gctl::callink_magnetic_para_earth_sph(...)");
|
||||
}
|
||||
|
||||
point3dc v1, v2, v3, ne, nf, mag_z, c;
|
||||
point3ds s;
|
||||
|
||||
double I = inclina_deg*M_PI/180.0;
|
||||
double A = declina_deg*M_PI/180.0;
|
||||
|
||||
if (mag_vec != nullptr) mag_vec->resize(in_cone.size());
|
||||
out_para.resize(in_cone.size());
|
||||
for (int i = 0; i < in_cone.size(); ++i)
|
||||
{
|
||||
if (in_cone[i].vert[0] == nullptr || in_cone[i].vert[1] == nullptr ||
|
||||
in_cone[i].vert[2] == nullptr || in_cone[i].vert[3] == nullptr)
|
||||
{
|
||||
throw runtime_error("Invalid vertex pointer. From callink_magnetic_para_earth_sph(...)");
|
||||
}
|
||||
|
||||
// magnetization vector at the local coordinates
|
||||
// note here the postive direction of the inclination angle is downward
|
||||
// note here the postive direction of the declination angle is clockwise
|
||||
mag_z.x = cos(I)*sin(A);
|
||||
mag_z.y = cos(I)*cos(A);
|
||||
mag_z.z = -1.0*sin(I);
|
||||
|
||||
c = 1.0/3.0*(*in_cone[i].vert[0] + *in_cone[i].vert[1] + *in_cone[i].vert[2]);
|
||||
s = c.c2s();
|
||||
// magnetic susceptibility is taken as one here
|
||||
// rotate the local coordinate system to the regular status
|
||||
mag_z = field_tense * mag_z.rotate((90.0 - s.lat)*M_PI/180.0, 0.0, (90.0 + s.lon)*M_PI/180.0).normal();
|
||||
if (mag_vec != nullptr) mag_vec->at(i) = mag_z;
|
||||
|
||||
for (int f = 0; f < 4; ++f)
|
||||
{
|
||||
v1 = *in_cone[i].fget(f, 1) - *in_cone[i].fget(f, 0);
|
||||
v2 = *in_cone[i].fget(f, 2) - *in_cone[i].fget(f, 0);
|
||||
nf = cross(v1, v2).normal();
|
||||
out_para[i].nf[f] = nf;
|
||||
|
||||
for (int e = 0; e < 3; ++e)
|
||||
{
|
||||
v3 = *in_cone[i].fget(f, (e+1)%3) - *in_cone[i].fget(f, e);
|
||||
ne = cross(v3, nf).normal();
|
||||
out_para[i].ne[e+f*3] = ne;
|
||||
out_para[i].te[e+f*3] = cross(nf, ne).normal();
|
||||
}
|
||||
|
||||
out_para[i].mag_amp[f] = dot(mag_z, nf);
|
||||
}
|
||||
|
||||
in_cone[i].att = out_para.get(i);
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
gctl::point3dc gctl::magkernel_single(const magcone_ren17 &a_ele, const point3ds &a_op, tensor *R_ptr)
|
||||
{
|
||||
double Rij_minus, Rij_plus, Sij_plus, Sij_minus, Rij0, mij0, wi0;
|
||||
double beta, Aij, sig, absw;
|
||||
gctl::point3dc oi, k1, part1, part2;
|
||||
|
||||
// get attribute pointer
|
||||
if (a_ele.att == nullptr) throw gctl::runtime_error("[gctl_potential::magcone_ren17] Magnetization parameter not set.");
|
||||
gctl::magcone_para_ren17 *mpara = a_ele.att;
|
||||
|
||||
gctl::tensor R;
|
||||
if (R_ptr != nullptr) R = *R_ptr;
|
||||
else R = transform_matrix(a_op);
|
||||
|
||||
// get the observation site in the local coordinates
|
||||
gctl::point3dc site = a_op.s2c();
|
||||
gctl::point3dc out_grad(0.0, 0.0, 0.0);
|
||||
for (int f = 0; f < 4; ++f)
|
||||
{
|
||||
k1.set(0.0, 0.0, 0.0);
|
||||
for (int j = 0; j < 3; ++j)
|
||||
{
|
||||
Rij_minus = (site - *a_ele.fget(f, j)).module();
|
||||
Rij_plus = (site - *a_ele.fget(f, (j+1)%3)).module();
|
||||
|
||||
if (j == 0)
|
||||
{
|
||||
wi0 = gctl::dot(site - *a_ele.fget(f, j), mpara->nf[f]);
|
||||
sig = gctl::sign(wi0);
|
||||
absw = std::abs(wi0);
|
||||
}
|
||||
|
||||
oi = site - wi0*mpara->nf[f];
|
||||
Sij_minus = gctl::dot(*a_ele.fget(f, j) - oi, mpara->te[3*f+j]);
|
||||
Sij_plus = gctl::dot(*a_ele.fget(f, (j+1)%3) - oi, mpara->te[3*f+j]);
|
||||
mij0 = gctl::dot(*a_ele.fget(f, j) - oi, mpara->ne[3*f+j]);
|
||||
Rij0 = std::sqrt(wi0*wi0 + mij0*mij0);
|
||||
|
||||
part2.set(0.0, 0.0, 0.0);
|
||||
if (absw > mag_tolerance)
|
||||
{
|
||||
beta = atan((mij0*Sij_plus)/(Rij0*Rij0 + absw*Rij_plus))
|
||||
- atan((mij0*Sij_minus)/(Rij0*Rij0 + absw*Rij_minus));
|
||||
|
||||
part2 = sig*beta*mpara->nf[f];
|
||||
}
|
||||
|
||||
if (std::abs(Rij0) > mag_tolerance)
|
||||
{
|
||||
Aij = std::log((long double)(Rij_plus+Sij_plus)) - std::log((long double)(Rij_minus+Sij_minus));
|
||||
}
|
||||
else if (Sij_plus > 0.0 && Sij_minus > 0.0)
|
||||
{
|
||||
Aij = std::log(Sij_plus) - std::log(Sij_minus);
|
||||
}
|
||||
else if (Sij_plus < 0.0 && Sij_minus < 0.0)
|
||||
{
|
||||
Aij = std::log(-1.0*Sij_minus) - std::log(-1.0*Sij_plus);
|
||||
}
|
||||
else
|
||||
{
|
||||
throw gctl::runtime_error("[gctl_potential::magcone_ren17] Observation site on edge.");
|
||||
}
|
||||
|
||||
part1 = Aij * mpara->ne[3*f+j];
|
||||
k1 = k1 - (part1 + part2);
|
||||
}
|
||||
|
||||
// nT -> T 1e-9 / 4*pi*e-7 = 1/(400*pi)*100 = 1/(4*pi) -> A/m
|
||||
out_grad = out_grad - 1.0/(4.0*GCTL_Pi)*k1*mpara->mag_amp[f];
|
||||
}
|
||||
|
||||
out_grad = R*out_grad;
|
||||
out_grad.x *= -1.0;
|
||||
out_grad.y *= -1.0;
|
||||
return out_grad;
|
||||
}
|
||||
|
||||
gctl::tensor gctl::magkernel_single_tensor(const magcone_ren17 &a_ele, const point3ds &a_op, tensor *R_ptr)
|
||||
{
|
||||
double Rij_minus, Rij_plus, Sij_plus, Sij_minus, Rij0, mij0, wi0;
|
||||
double beta, sig, absw;
|
||||
double factor_n_mij, factor_tij;
|
||||
gctl::point3dc oi, grad_Aij;
|
||||
gctl::point3dc grad_Rij_plus, grad_Rij_minus, grad_Sij_plus, grad_Sij_minus;
|
||||
gctl::point3dc grad_mij0, grad_Rij0, grad_abs_wi0;
|
||||
gctl::point3dc grad_a_plus, grad_b_plus, grad_a_minus, grad_b_minus;
|
||||
gctl::point3dc grad_betaij_plus, grad_betaij_minus, grad_betaij;
|
||||
double a_plus, b_plus, a_minus, b_minus;
|
||||
double k3;
|
||||
gctl::tensor tmp_k, k2;
|
||||
|
||||
// get attribute pointer
|
||||
if (a_ele.att == nullptr) throw gctl::runtime_error("[gctl_potential::magcone_ren17] Magnetization parameter not set.");
|
||||
gctl::magcone_para_ren17 *mpara = a_ele.att;
|
||||
|
||||
gctl::tensor R, R_T;
|
||||
if (R_ptr != nullptr) R = *R_ptr;
|
||||
else R = transform_matrix(a_op);
|
||||
R_T = R.transpose();
|
||||
|
||||
// get the observation site in the local coordinates
|
||||
gctl::point3dc site = a_op.s2c();
|
||||
gctl::tensor out_tensor(0.0);
|
||||
for (int f = 0; f < 4; ++f)
|
||||
{
|
||||
k2.set(0.0);
|
||||
k3 = 0.0;
|
||||
for (int j = 0; j < 3; ++j)
|
||||
{
|
||||
Rij_minus = (site - *a_ele.fget(f, j)).module();
|
||||
Rij_plus = (site - *a_ele.fget(f, (j+1)%3)).module();
|
||||
|
||||
if (j == 0)
|
||||
{
|
||||
wi0 = gctl::dot(site - *a_ele.fget(f, j), mpara->nf[f]);
|
||||
sig = gctl::sign(wi0);
|
||||
absw = std::abs(wi0);
|
||||
}
|
||||
|
||||
oi = site - wi0*mpara->nf[f];
|
||||
Sij_minus = gctl::dot(*a_ele.fget(f, j) - oi, mpara->te[3*f+j]);
|
||||
Sij_plus = gctl::dot(*a_ele.fget(f, (j+1)%3) - oi, mpara->te[3*f+j]);
|
||||
mij0 = gctl::dot(*a_ele.fget(f, j) - oi, mpara->ne[3*f+j]);
|
||||
Rij0 = std::sqrt(wi0*wi0 + mij0*mij0);
|
||||
|
||||
if (std::abs(Rij0) > mag_tolerance)
|
||||
{
|
||||
factor_n_mij = -1.0*(Sij_plus/(Rij0*Rij0*Rij_plus) - Sij_minus/(Rij0*Rij0*Rij_minus));
|
||||
factor_tij = -1.0/Rij_plus + 1.0/Rij_minus;
|
||||
grad_Aij = (wi0*factor_n_mij)*mpara->nf[f] + factor_tij*mpara->te[3*f+j]
|
||||
- (mij0*factor_n_mij)*mpara->ne[3*f+j];
|
||||
}
|
||||
else
|
||||
{
|
||||
factor_tij = -1.0/Rij_plus + 1.0/Rij_minus;
|
||||
grad_Aij = factor_tij*mpara->te[3*f+j];
|
||||
}
|
||||
|
||||
//tmp_k = gctl::kron(grad_Aij, ne_[e][3*f+j]);
|
||||
tmp_k = gctl::kron(mpara->ne[3*f+j], grad_Aij);
|
||||
k2 = k2 - tmp_k;
|
||||
|
||||
if (absw > mag_tolerance)
|
||||
{
|
||||
grad_Rij_plus = (1.0/Rij_plus)*(site - *a_ele.fget(f, (j+1)%3));
|
||||
grad_Rij_minus = (1.0/Rij_minus)*(site - *a_ele.fget(f, j));
|
||||
grad_Sij_plus = -1.0*mpara->te[3*f+j];
|
||||
grad_Sij_minus = -1.0*mpara->te[3*f+j];
|
||||
grad_mij0 = -1.0*mpara->ne[3*f+j];
|
||||
grad_Rij0 = (1.0/Rij0)*(wi0*mpara->nf[f] - mij0*mpara->ne[3*f+j]);
|
||||
grad_abs_wi0 = sig*mpara->nf[f];
|
||||
a_plus = Rij0*Rij0 + absw*Rij_plus;
|
||||
b_plus = mij0*Sij_plus;
|
||||
grad_a_plus = (2.0*Rij0)*grad_Rij0 + Rij_plus*grad_abs_wi0 + absw*grad_Rij_plus;
|
||||
grad_b_plus = Sij_plus*grad_mij0 + mij0*grad_Sij_plus;
|
||||
a_minus = Rij0*Rij0 + absw*Rij_minus;
|
||||
b_minus = mij0*Sij_minus;
|
||||
grad_a_minus = (2.0*Rij0)*grad_Rij0 + Rij_minus*grad_abs_wi0 + absw*grad_Rij_minus;
|
||||
grad_b_minus = Sij_minus*grad_mij0 + mij0*grad_Sij_minus;
|
||||
grad_betaij_plus = (1.0/(a_plus*a_plus + b_plus*b_plus))*(a_plus*grad_b_plus - b_plus*grad_a_plus);
|
||||
grad_betaij_minus = (1.0/(a_minus*a_minus + b_minus*b_minus))*(a_minus*grad_b_minus - b_minus*grad_a_minus);
|
||||
|
||||
grad_betaij = grad_betaij_plus - grad_betaij_minus;
|
||||
|
||||
tmp_k = gctl::kron(grad_betaij, mpara->nf[f]);
|
||||
k2 = k2 - sig*tmp_k;
|
||||
}
|
||||
else
|
||||
{
|
||||
if (std::abs(mij0) > mag_tolerance)
|
||||
{
|
||||
k3 += (-1.0/mij0)*(Sij_plus/Rij_plus - Sij_minus/Rij_minus);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (k3 != 0.0)
|
||||
{
|
||||
tmp_k = gctl::kron(mpara->nf[f], mpara->nf[f]);
|
||||
k2 = k2 - k3*tmp_k;
|
||||
}
|
||||
|
||||
// nT -> T 1e-9 / 4*pi*e-7 = 1/(400*pi)*100 = 1/(4*pi) -> A/m
|
||||
out_tensor = out_tensor - 1.0/(4.0*GCTL_Pi)*k2*mpara->mag_amp[f];
|
||||
}
|
||||
|
||||
out_tensor = R*out_tensor*R_T;
|
||||
out_tensor[0][0] *= -1.0;
|
||||
out_tensor[0][1] *= -1.0;
|
||||
out_tensor[0][2] *= -1.0;
|
||||
out_tensor[1][0] *= -1.0;
|
||||
out_tensor[2][0] *= -1.0;
|
||||
return out_tensor;
|
||||
}
|
||||
|
||||
void gctl::magkernel(matrix<double> &kernel, const array<magcone_ren17> &top_ele,
|
||||
const array<magcone_ren17> &btm_ele, const array<point3ds> &obsp,
|
||||
magnetic_field_type_e comp_type, verbose_type_e verbose)
|
||||
{
|
||||
if (comp_type != Bx && comp_type != By && comp_type != Bz)
|
||||
{
|
||||
throw std::runtime_error("[gctl_potential::magcone_ren17] Invalid magnetic componment type. Must be Bx, By or Bz.");
|
||||
return;
|
||||
}
|
||||
|
||||
if (top_ele.size() != btm_ele.size())
|
||||
{
|
||||
throw std::runtime_error("[gctl_potential::magcone_ren17] Elements' size don't equal.");
|
||||
return;
|
||||
}
|
||||
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = top_ele.size();
|
||||
kernel.resize(o_size, e_size);
|
||||
|
||||
array<tensor> Rs(o_size);
|
||||
#pragma omp parallel for private (i) schedule(guided)
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
Rs[i] = transform_matrix(obsp[i]);
|
||||
}
|
||||
|
||||
point3dc mag_b;
|
||||
gctl::progress_bar bar(o_size*2, "magkernel_componments");
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(i);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);
|
||||
|
||||
#pragma omp parallel for private (j, mag_b) shared(i) schedule(guided)
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
mag_b = magkernel_single(top_ele[j], obsp[i], Rs.get(i));
|
||||
if (comp_type == Bx) kernel[i][j] = mag_b.y;
|
||||
if (comp_type == By) kernel[i][j] = mag_b.z;
|
||||
if (comp_type == Bz) kernel[i][j] = mag_b.x;
|
||||
}
|
||||
}
|
||||
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(i + o_size);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(i + o_size);
|
||||
|
||||
#pragma omp parallel for private (j, mag_b) shared(i) schedule(guided)
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
mag_b = magkernel_single(btm_ele[j], obsp[i], Rs.get(i));
|
||||
if (comp_type == Bx) kernel[i][j] -= mag_b.y;
|
||||
if (comp_type == By) kernel[i][j] -= mag_b.z;
|
||||
if (comp_type == Bz) kernel[i][j] -= mag_b.x;
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::magkernel(spmat<double> &kernel, const array<magcone_ren17> &top_ele,
|
||||
const array<magcone_ren17> &btm_ele, const array<point3ds> &obsp,
|
||||
double cut_angle, magnetic_field_type_e comp_type, verbose_type_e verbose)
|
||||
{
|
||||
if (comp_type != Bx && comp_type != By && comp_type != Bz)
|
||||
{
|
||||
throw std::runtime_error("[gctl_potential::magcone_ren17] Invalid magnetic componment type. Must be Bx, By or Bz.");
|
||||
return;
|
||||
}
|
||||
|
||||
if (top_ele.size() != btm_ele.size())
|
||||
{
|
||||
throw std::runtime_error("[gctl_potential::magcone_ren17] Elements' size don't equal.");
|
||||
return;
|
||||
}
|
||||
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = top_ele.size();
|
||||
|
||||
array<tensor> Rs(o_size);
|
||||
#pragma omp parallel for private (i) schedule(guided)
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
Rs[i] = transform_matrix(obsp[i]);
|
||||
}
|
||||
|
||||
point3dc cen, mag_b;
|
||||
mat_node<double> tmp_plt;
|
||||
std::map<int, int> tri_idx;
|
||||
std::vector<mat_node<double>> triplts;
|
||||
|
||||
gctl::progress_bar bar2(o_size, "initializing triplts");
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar2.progressed(i);
|
||||
else if (verbose == gctl::ShortMsg) bar2.progressed_simple(i);
|
||||
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
cen = 1.0/3.0*(*top_ele[j].vert[0] + *top_ele[j].vert[1] + *top_ele[j].vert[2]);
|
||||
if (geometry3d::angle(obsp[i].s2c(), cen) < cut_angle*GCTL_Pi/180.0)
|
||||
{
|
||||
tmp_plt.r_id = i;
|
||||
tmp_plt.c_id = j;
|
||||
tmp_plt.val = 1.0;
|
||||
|
||||
tri_idx[j + i*e_size] = triplts.size();
|
||||
triplts.push_back(tmp_plt);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
gctl::progress_bar bar(o_size*2, "magkernel_componments");
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(i);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);
|
||||
|
||||
#pragma omp parallel for private (j, mag_b, cen) shared(i, cut_angle) schedule(guided)
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
cen = 1.0/3.0*(*top_ele[j].vert[0] + *top_ele[j].vert[1] + *top_ele[j].vert[2]);
|
||||
if (geometry3d::angle(obsp[i].s2c(), cen) < cut_angle*GCTL_Pi/180.0)
|
||||
{
|
||||
mag_b = magkernel_single(top_ele[j], obsp[i], Rs.get(i));
|
||||
if (comp_type == Bx) triplts[tri_idx[j + i*e_size]].val = mag_b.y;
|
||||
if (comp_type == By) triplts[tri_idx[j + i*e_size]].val = mag_b.z;
|
||||
if (comp_type == Bz) triplts[tri_idx[j + i*e_size]].val = mag_b.x;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(i + o_size);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(i + o_size);
|
||||
|
||||
#pragma omp parallel for private (j, mag_b, cen) shared(i, cut_angle) schedule(guided)
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
cen = 1.0/3.0*(*btm_ele[j].vert[0] + *btm_ele[j].vert[1] + *btm_ele[j].vert[2]);
|
||||
if (geometry3d::angle(obsp[i].s2c(), cen) < cut_angle*GCTL_Pi/180.0)
|
||||
{
|
||||
mag_b = magkernel_single(btm_ele[j], obsp[i], Rs.get(i));
|
||||
if (comp_type == Bx) triplts[tri_idx[j + i*e_size]].val -= mag_b.y;
|
||||
if (comp_type == By) triplts[tri_idx[j + i*e_size]].val -= mag_b.z;
|
||||
if (comp_type == Bz) triplts[tri_idx[j + i*e_size]].val -= mag_b.x;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
kernel.malloc(o_size, e_size, 0.0);
|
||||
kernel.set_triplts(triplts);
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::magobser(array<point3dc> &out_obs, const array<magcone_ren17> &top_ele,
|
||||
const array<magcone_ren17> &btm_ele, const array<point3ds> &obsp,
|
||||
const array<double> &sus, verbose_type_e verbose)
|
||||
{
|
||||
if (top_ele.size() != btm_ele.size())
|
||||
{
|
||||
throw std::runtime_error("[gctl_potential::magcone_ren17] Elements' size don't equal.");
|
||||
return;
|
||||
}
|
||||
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = top_ele.size();
|
||||
|
||||
out_obs.resize(o_size, point3dc(0.0, 0.0, 0.0));
|
||||
array<tensor> Rs(o_size);
|
||||
#pragma omp parallel for private (i) schedule(guided)
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
Rs[i] = transform_matrix(obsp[i]);
|
||||
}
|
||||
|
||||
gctl::progress_bar bar(e_size*2, "magobser_componments");
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(j);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
|
||||
|
||||
#pragma omp parallel for private (i) shared (j) schedule(guided)
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
out_obs[i] = out_obs[i] + sus[j] * magkernel_single(top_ele[j], obsp[i], Rs.get(i));
|
||||
}
|
||||
}
|
||||
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(j + e_size);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(j + e_size);
|
||||
|
||||
#pragma omp parallel for private (i) shared (j) schedule(guided)
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
out_obs[i] = out_obs[i] - sus[j] * magkernel_single(btm_ele[j], obsp[i], Rs.get(i));
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
void gctl::magobser(array<tensor> &out_obs, const array<magcone_ren17> &top_ele,
|
||||
const array<magcone_ren17> &btm_ele, const array<point3ds> &obsp,
|
||||
const array<double> &sus, verbose_type_e verbose)
|
||||
{
|
||||
if (top_ele.size() != btm_ele.size())
|
||||
{
|
||||
throw std::runtime_error("[gctl_potential::magcone_ren17] Elements' size don't equal.");
|
||||
return;
|
||||
}
|
||||
|
||||
int i, j;
|
||||
int o_size = obsp.size();
|
||||
int e_size = top_ele.size();
|
||||
|
||||
out_obs.resize(o_size, tensor(0.0));
|
||||
array<tensor> Rs(o_size);
|
||||
#pragma omp parallel for private (i) schedule(guided)
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
Rs[i] = transform_matrix(obsp[i]);
|
||||
}
|
||||
|
||||
gctl::progress_bar bar(e_size*2, "magobser_componments");
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(j);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
|
||||
|
||||
#pragma omp parallel for private (i) shared (j) schedule(guided)
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
out_obs[i] = out_obs[i] + sus[j] * magkernel_single_tensor(top_ele[j], obsp[i], Rs.get(i));
|
||||
}
|
||||
}
|
||||
|
||||
for (j = 0; j < e_size; j++)
|
||||
{
|
||||
if (verbose == gctl::FullMsg) bar.progressed(j + e_size);
|
||||
else if (verbose == gctl::ShortMsg) bar.progressed_simple(j + e_size);
|
||||
|
||||
#pragma omp parallel for private (i) shared (j) schedule(guided)
|
||||
for (i = 0; i < o_size; i++)
|
||||
{
|
||||
out_obs[i] = out_obs[i] - sus[j] * magkernel_single_tensor(btm_ele[j], obsp[i], Rs.get(i));
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
142
lib/potential/mkernel_tricone_Ren2017.h
Normal file
142
lib/potential/mkernel_tricone_Ren2017.h
Normal file
@ -0,0 +1,142 @@
|
||||
/********************************************************
|
||||
* ██████╗ ██████╗████████╗██╗
|
||||
* ██╔════╝ ██╔════╝╚══██╔══╝██║
|
||||
* ██║ ███╗██║ ██║ ██║
|
||||
* ██║ ██║██║ ██║ ██║
|
||||
* ╚██████╔╝╚██████╗ ██║ ███████╗
|
||||
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
|
||||
* Geophysical Computational Tools & Library (GCTL)
|
||||
*
|
||||
* Copyright (c) 2022 Yi Zhang (yizhang-geo@zju.edu.cn)
|
||||
*
|
||||
* GCTL is distributed under a dual licensing scheme. You can redistribute
|
||||
* it and/or modify it under the terms of the GNU Lesser General Public
|
||||
* License as published by the Free Software Foundation, either version 2
|
||||
* of the License, or (at your option) any later version. You should have
|
||||
* received a copy of the GNU Lesser General Public License along with this
|
||||
* program. If not, see <http://www.gnu.org/licenses/>.
|
||||
*
|
||||
* If the terms and conditions of the LGPL v.2. would prevent you from using
|
||||
* the GCTL, please consider the option to obtain a commercial license for a
|
||||
* fee. These licenses are offered by the GCTL's original author. As a rule,
|
||||
* licenses are provided "as-is", unlimited in time for a one time fee. Please
|
||||
* send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget
|
||||
* to include some description of your company and the realm of its activities.
|
||||
* Also add information on how to contact you by electronic and paper mail.
|
||||
******************************************************/
|
||||
|
||||
#ifndef _GCTL_MAG_KERNEL_TRICONE_REN2017_H
|
||||
#define _GCTL_MAG_KERNEL_TRICONE_REN2017_H
|
||||
|
||||
#include "gm_data.h"
|
||||
|
||||
namespace gctl
|
||||
{
|
||||
struct magcone_para_ren17
|
||||
{
|
||||
double mag_amp[4]; // 四面体四个面的磁化强度
|
||||
point3dc nf[4]; // 四面体面外法线矢量
|
||||
point3dc ne[12]; // 四面体边外法线矢量
|
||||
point3dc te[12]; // 四面体边切线矢量
|
||||
};
|
||||
|
||||
typedef type_tricone<magcone_para_ren17> magcone_ren17; ///< 带magcone_para_ren17属性的三角锥结构体
|
||||
|
||||
/**
|
||||
* @brief Set the tolerance for calculating the magnetic parameters of tricone elements.
|
||||
*
|
||||
* @param tol Tolerance
|
||||
*/
|
||||
void set_magcone_ren17_tolerance(double tol);
|
||||
|
||||
/**
|
||||
* @brief Calculate the magnetic parameters of given tricone elements.
|
||||
*
|
||||
* @param in_tet Input and output elements
|
||||
* @param out_para Output parameters
|
||||
* @param mag_B magnetization vecrtors
|
||||
*/
|
||||
void callink_magnetic_para(array<magcone_ren17> &in_cone, array<magcone_para_ren17> &out_para, const array<point3dc> &mag_B);
|
||||
|
||||
/**
|
||||
* @brief Calculate the magnetic parameters of given tricone elements wrt. the spherical coordinates.
|
||||
*
|
||||
* @note The value of magnetic susceptibility is taken as one here. This is usefull for calculating
|
||||
* kernel matrix of the magnetic anomalies.
|
||||
*
|
||||
* @param in_tet Input elements
|
||||
* @param out_para Output parameters
|
||||
* @param inclina_deg inclination angle of the magnetization vector wrt. the local Cartesian coordinates at every tricone elements
|
||||
* @param declina_deg declination angle of the magnetization vector wrt. the local Cartesian coordinates at every tricone elements
|
||||
* @param mag_vec Output magnetization vectors (This is useful for data visualization)
|
||||
* @param field_tense Tense of the Earth's magnetic field
|
||||
*/
|
||||
void callink_magnetic_para_earth_sph(array<magcone_ren17> &in_tet, array<magcone_para_ren17> &out_para,
|
||||
double inclina_deg, double declina_deg, array<point3dc> *mag_vec = nullptr, double field_tense = GCTL_T0);
|
||||
|
||||
/**
|
||||
* @brief Calculate the magnetic field vector at a given observation point.
|
||||
*
|
||||
* @param a_ele Input element
|
||||
* @param a_op Input observation point
|
||||
* @param R_ptr Output rotation matrix
|
||||
* @return point3dc Magnetic field vector
|
||||
*/
|
||||
point3dc magkernel_single(const magcone_ren17 &a_ele, const point3ds &a_op, tensor *R_ptr = nullptr);
|
||||
|
||||
/**
|
||||
* @brief Calculate the magnetic field vector at a given observation point.
|
||||
*
|
||||
* @param a_ele Input element
|
||||
* @param a_op Input observation point
|
||||
* @param R_ptr Output rotation matrix
|
||||
* @return point3dc Magnetic field vector
|
||||
*/
|
||||
tensor magkernel_single_tensor(const magcone_ren17 &a_ele, const point3ds &a_op, tensor *R_ptr = nullptr);
|
||||
|
||||
/**
|
||||
* @brief Calculate the magnetic field vector at a given observation point.
|
||||
*
|
||||
* @param a_ele Input element
|
||||
* @param a_op Input observation point
|
||||
* @param R_ptr Output rotation matrix
|
||||
* @return point3dc Magnetic field vector
|
||||
*/
|
||||
void magkernel(matrix<double> &kernel, const array<magcone_ren17> &top_ele, const array<magcone_ren17> &btm_ele,
|
||||
const array<point3ds> &obsp, magnetic_field_type_e comp_type = Bz, verbose_type_e verbose = FullMsg);
|
||||
|
||||
/**
|
||||
* @brief Calculate the magnetic field vector at a given observation point.
|
||||
*
|
||||
* @param a_ele Input element
|
||||
* @param a_op Input observation point
|
||||
* @param R_ptr Output rotation matrix
|
||||
* @return point3dc Magnetic field vector
|
||||
*/
|
||||
void magkernel(spmat<double> &kernel, const array<magcone_ren17> &top_ele, const array<magcone_ren17> &btm_ele,
|
||||
const array<point3ds> &obsp, double cut_angle, magnetic_field_type_e comp_type = Bz, verbose_type_e verbose = FullMsg);
|
||||
|
||||
/**
|
||||
* @brief Calculate the magnetic field vector at a given observation point.
|
||||
*
|
||||
* @param a_ele Input element
|
||||
* @param a_op Input observation point
|
||||
* @param R_ptr Output rotation matrix
|
||||
* @return point3dc Magnetic field vector
|
||||
*/
|
||||
void magobser(array<point3dc> &out_obs, const array<magcone_ren17> &top_ele, const array<magcone_ren17> &btm_ele,
|
||||
const array<point3ds> &obsp, const array<double> &sus, verbose_type_e verbose = FullMsg);
|
||||
|
||||
/**
|
||||
* @brief Calculate the magnetic field vector at a given observation point.
|
||||
*
|
||||
* @param a_ele Input element
|
||||
* @param a_op Input observation point
|
||||
* @param R_ptr Output rotation matrix
|
||||
* @return point3dc Magnetic field vector
|
||||
*/
|
||||
void magobser(array<tensor> &out_obs, const array<magcone_ren17> &top_ele, const array<magcone_ren17> &btm_ele,
|
||||
const array<point3ds> &obsp, const array<double> &sus, verbose_type_e verbose = FullMsg);
|
||||
}
|
||||
|
||||
#endif // _GCTL_MAG_KERNEL_TRICONE_REN2017_H
|
||||
@ -372,8 +372,7 @@ void save_gmsh(const std::vector<std::string> &cmd_units)
|
||||
// save gmsh <file>
|
||||
if (cmd_units.size() < 3) throw std::runtime_error("save: insufficient parameters.");
|
||||
|
||||
rg.save_gmsh(cmd_units[2], NodeData, OverWrite, NotPacked);
|
||||
rg.save_gmsh(cmd_units[2], ElemData, Append, NotPacked);
|
||||
rg.save_gmsh(cmd_units[2], OverWrite, NotPacked);
|
||||
return;
|
||||
}
|
||||
|
||||
@ -472,7 +471,7 @@ void data_cloud(const std::vector<std::string> &cmd_units)
|
||||
array<point2dc> posi_arr(posix_vec.size());
|
||||
array<double> posi_val;
|
||||
|
||||
posi_val.import_vector(data_vec);
|
||||
posi_val.input(data_vec);
|
||||
for (size_t i = 0; i < posi_arr.size(); i++)
|
||||
{
|
||||
posi_arr[i].x = posix_vec[i];
|
||||
@ -622,16 +621,16 @@ void data_output(const std::vector<std::string> &cmd_units)
|
||||
{
|
||||
for (size_t i = 1; i < copy_str.size(); i++)
|
||||
{
|
||||
data_ptr = rg.get_data(copy_str[i]);
|
||||
data_ptr->set_output(true);
|
||||
data_ptr = rg.get_data_ptr(copy_str[i]);
|
||||
data_ptr->output_ok_ = true;
|
||||
}
|
||||
}
|
||||
else if (copy_str[0] == "disable")
|
||||
{
|
||||
for (size_t i = 1; i < copy_str.size(); i++)
|
||||
{
|
||||
data_ptr = rg.get_data(copy_str[i]);
|
||||
data_ptr->set_output(false);
|
||||
data_ptr = rg.get_data_ptr(copy_str[i]);
|
||||
data_ptr->output_ok_ = false;
|
||||
}
|
||||
}
|
||||
else throw std::runtime_error("data-output: invalid operation type.");
|
||||
@ -649,7 +648,8 @@ void data_rename(const std::vector<std::string> &cmd_units)
|
||||
copy_str[i - 1] = cmd_units[i];
|
||||
}
|
||||
|
||||
rg.rename_data(copy_str[0], copy_str[1]);
|
||||
meshdata &data = rg.get_data(copy_str[0]);
|
||||
data.name_ = copy_str[1];
|
||||
return;
|
||||
}
|
||||
|
||||
@ -692,7 +692,7 @@ void get_stats(const std::vector<std::string> &cmd_units)
|
||||
std::vector<double> stats;
|
||||
for (size_t i = 0; i < copy_str.size(); i++)
|
||||
{
|
||||
data_ptr = rg.get_data(copy_str[i]);
|
||||
data_ptr = rg.get_data_ptr(copy_str[i]);
|
||||
data_ptr->show_info();
|
||||
data_ptr->show_stats();
|
||||
}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user