Skip to content

GitLab

Commit 66a48942 authored by Milad Bader's avatar Milad Bader
Browse files

Delta function discretization with 3 moment conditions added

parent 4b4a9e0f
master
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 458 additions and 42 deletions
+458 -42
src/cpp/main/elastic_we_modeling.cpp
+ 4
- 4
View file @ 66a48942
......@@ -138,7 +138,7 @@ int main(int argc, char **argv) {
MB::readParam(argc, argv, "regrid", regrid,false);
MB::readParam(argc, argv, "resamp_type", resamp_type,"linear");
MB::readParam(argc, argv, "sinc_half_length", sinc_half_length,40);
MB::readParam(argc, argv, "injector_type", injector_type,"nearest_neighbor");
MB::readParam(argc, argv, "injector_type", injector_type,"delta_m3");
MB::readParam(argc, argv, "wavefield", wavefield, "none");
MB::readParam(argc, argv, "info", info, false);
MB::readParam(argc, argv, "ispc", ispc, false);
......@@ -171,11 +171,11 @@ int main(int argc, char **argv) {
bool moment_tensor = false;
std::string rcv_injector_type=injector_type;
if ((srcxz_file != "none") || (injector_type == "delta_derivative_m3") || (injector_type == "delta_derivative_m5"))
if ((injector_type == "delta_derivative_m3") || (injector_type == "delta_derivative_m5"))
moment_tensor=true;
if (moment_tensor == true) {
if ((injector_type == "delta_derivative_m3") || (injector_type == "delta_m2")) {rcv_injector_type = "delta_m2"; injector_type="delta_derivative_m3";}
if ((injector_type == "delta_derivative_m5") || (injector_type == "delta_m4")) {rcv_injector_type = "delta_m4"; injector_type="delta_derivative_m5";}
if (injector_type == "delta_derivative_m3") {rcv_injector_type = "delta_m3"; injector_type="delta_derivative_m3";}
if (injector_type == "delta_derivative_m5") {rcv_injector_type = "delta_m4"; injector_type="delta_derivative_m5";}
}
......
src/cpp/scalar/injecter.cpp
+ 219
- 23
View file @ 66a48942
......@@ -284,6 +284,118 @@ void injecter_delta_m2::adjoint(bool add, std::shared_ptr<vec2DReg> mod, const s
}
void injecter_delta_m3::findIndex(const std::shared_ptr<hypercube> range){
data_t dx = range->getAxis(2).d;
data_t dz = range->getAxis(1).d;
int nx = range->getAxis(2).n;
int nz = range->getAxis(1).n;
int ox = range->getAxis(2).o;
int oz = range->getAxis(1).o;
for (int isr = 0; isr < _nsr; isr++){
assert( (_xloc[isr] >= ox) && (_xloc[isr] <= ox+(nx-1)*dx) && (_zloc[isr] >= oz) && (_zloc[isr] <= oz+(nz-1)*dz) );
// compute the index of the pivot
_xIndex[isr][1] = std::floor((_xloc[isr] - ox)/dx);
_zIndex[isr][1] = std::floor((_zloc[isr] - oz)/dz);
// check if the pivot falls at or near the boundary
// if so, shift it forward or backward one or two samples to make room for the other points xi-1 xi+1 and xi+2
if (_xIndex[isr][1] == 0) _xIndex[isr][1] = 1;
if (_xIndex[isr][1] == nx - 1) _xIndex[isr][1] = nx - 2;
if (_zIndex[isr][1] == 0) _zIndex[isr][1] = 1;
if (_zIndex[isr][1] == nz - 1) _zIndex[isr][1] = nz - 2;
_xIndex[isr][0] = _xIndex[isr][1] - 1;
_xIndex[isr][2] = _xIndex[isr][1] + 1;
_zIndex[isr][0] = _zIndex[isr][1] - 1;
_zIndex[isr][2] = _zIndex[isr][1] + 1;
}
}
void injecter_delta_m3::findWeight(const HxHz * H){
data_t dx = _range->getAxis(2).d;
data_t dz = _range->getAxis(1).d;
int ox = _range->getAxis(2).o;
int oz = _range->getAxis(1).o;
data_t alpha_x, alpha_z;
for (int isr = 0; isr < _nsr; isr++){
// compute alpha based on the pivot: x = xi + alpha * h ; xi is the pivot
alpha_x = (_xloc[isr] - dx * _xIndex[isr][1] - ox) / dx;
alpha_z = (_zloc[isr] - dz * _zIndex[isr][1] - oz) / dz;
// compute the weights based on the formulas given in the header file
_xWeight[isr][0] = (pow(alpha_x,2) - alpha_x) / (2 * H->getWeightx(_xIndex[isr][0]));
_xWeight[isr][1] = ( -pow(alpha_x,2) + 1) / (H->getWeightx(_xIndex[isr][1]));
_xWeight[isr][2] = (pow(alpha_x,2) + alpha_x) / (2 * H->getWeightx(_xIndex[isr][2]));
_zWeight[isr][0] = (pow(alpha_z,2) - alpha_z) / (2 * H->getWeightz(_zIndex[isr][0]));
_zWeight[isr][1] = (-pow(alpha_z,2) + 1) / (H->getWeightz(_zIndex[isr][1]));
_zWeight[isr][2] = (pow(alpha_z,2) + alpha_z) / (2 * H->getWeightz(_zIndex[isr][2]));
_xWeight2[isr][0] = (pow(alpha_x,2) - alpha_x) / 2;
_xWeight2[isr][1] = ( -pow(alpha_x,2) + 1);
_xWeight2[isr][2] = (pow(alpha_x,2) + alpha_x) / 2;
_zWeight2[isr][0] = (pow(alpha_z,2) - alpha_z) / 2;
_zWeight2[isr][1] = ( -pow(alpha_z,2) + 1);
_zWeight2[isr][2] = (pow(alpha_z,2) + alpha_z) / 2;
}
}
void injecter_delta_m3::forward(bool add, const std::shared_ptr<vec2DReg> mod, std::shared_ptr<vec2DReg> dat) const{
this->checkDomainRange(mod,dat);
if (add == false) dat->zero();
int it = this->getIt();
int nz_mod = mod->getHyper()->getAxis(1).n;
int nz_dat = dat->getHyper()->getAxis(1).n;
for (int isr = 0; isr < _nsr; isr++){
for (int ix=0; ix<3; ix++){
for (int iz=0; iz<3; iz++){
dat->getVals()[_xIndex[isr][ix] * nz_dat + _zIndex[isr][iz]] += _xWeight[isr][ix] * _zWeight[isr][iz] * mod->getVals()[isr * nz_mod + it];
}
}
}
}
void injecter_delta_m3::adjoint(bool add, std::shared_ptr<vec2DReg> mod, const std::shared_ptr<vec2DReg> dat) const{
this->checkDomainRange(mod,dat);
if (add == false) mod->zero();
int it = this->getIt();
int nz_mod = mod->getHyper()->getAxis(1).n;
int nz_dat = dat->getHyper()->getAxis(1).n;
for (int isr = 0; isr < _nsr; isr++){
for (int ix=0; ix<3; ix++){
for (int iz=0; iz<3; iz++){
mod->getVals()[isr * nz_mod + it] += _xWeight2[isr][ix] * _zWeight2[isr][iz] * dat->getVals()[_xIndex[isr][ix] * nz_dat + _zIndex[isr][iz]];
}
}
}
}
void injecter_delta_m4::findIndex(const std::shared_ptr<hypercube> range){
data_t dx = range->getAxis(2).d;
......@@ -422,25 +534,13 @@ void injecter_delta_derivative_m3::findIndexDeriv(const std::shared_ptr<hypercub
_xIndexDeriv[isr][1] = std::floor((_xloc[isr] - ox)/dx);
_zIndexDeriv[isr][1] = std::floor((_zloc[isr] - oz)/dz);
// check if the pivot falls at the right/bottom boundary
// if so, shift it back one sample to make room for the other point xi+1
// check if the pivot falls at the boundary
// if so, shift it forward or backward one sample to make room for the other points xi-1 and xi+1
if (_xIndexDeriv[isr][1] == nx-1){
_xIndexDeriv[isr][1] = nx-2;
}
if (_zIndexDeriv[isr][1] == nz-1){
_zIndexDeriv[isr][1] = nz-2;
}
// check if the pivot falls at the left/top boundary
// if so, shift it forward one sample to make room for the other point xi-1
if (_xIndexDeriv[isr][1] == 0){
_xIndexDeriv[isr][1] = 1;
}
if (_zIndexDeriv[isr][1] == 0){
_zIndexDeriv[isr][1] = 1;
}
if (_xIndexDeriv[isr][1] == 0) _xIndexDeriv[isr][1] = 1;
if (_xIndexDeriv[isr][1] == nx - 1) _xIndexDeriv[isr][1] = nx - 2;
if (_zIndexDeriv[isr][1] == 0) _zIndexDeriv[isr][1] = 1;
if (_zIndexDeriv[isr][1] == nz - 1) _zIndexDeriv[isr][1] = nz - 2;
_xIndexDeriv[isr][0] = _xIndexDeriv[isr][1] - 1;
_zIndexDeriv[isr][0] = _zIndexDeriv[isr][1] - 1;
......@@ -490,7 +590,7 @@ void injecter_delta_derivative_m3::forward(bool add, const std::shared_ptr<vec2D
for (int isr = 0; isr < _nsr; isr++){
pdat[_xIndexDeriv[isr][0] * nz_dat + _zIndex[isr][0]] += _xWeightDeriv[isr][0] * _zWeight[isr][0] * pmod[isr * nz_mod + it];
/* pdat[_xIndexDeriv[isr][0] * nz_dat + _zIndex[isr][0]] += _xWeightDeriv[isr][0] * _zWeight[isr][0] * pmod[isr * nz_mod + it];
pdat[_xIndexDeriv[isr][1] * nz_dat + _zIndex[isr][0]] += _xWeightDeriv[isr][1] * _zWeight[isr][0] * pmod[isr * nz_mod + it];
pdat[_xIndexDeriv[isr][2] * nz_dat + _zIndex[isr][0]] += _xWeightDeriv[isr][2] * _zWeight[isr][0] * pmod[isr * nz_mod + it];
pdat[_xIndexDeriv[isr][0] * nz_dat + _zIndex[isr][1]] += _xWeightDeriv[isr][0] * _zWeight[isr][1] * pmod[isr * nz_mod + it];
......@@ -503,7 +603,13 @@ void injecter_delta_derivative_m3::forward(bool add, const std::shared_ptr<vec2D
pdat[_xIndex[isr][1] * nz_dat + _zIndexDeriv[isr][0]] += _xWeight[isr][1] * _zWeightDeriv[isr][0] * pmod[(_nsr+isr) * nz_mod + it];
pdat[_xIndex[isr][1] * nz_dat + _zIndexDeriv[isr][1]] += _xWeight[isr][1] * _zWeightDeriv[isr][1] * pmod[(_nsr+isr) * nz_mod + it];
pdat[_xIndex[isr][1] * nz_dat + _zIndexDeriv[isr][2]] += _xWeight[isr][1] * _zWeightDeriv[isr][2] * pmod[(_nsr+isr) * nz_mod + it];
*/
for (int j=0; j<3; j++){
for (int i=0; i<3; i++){
pdat[_xIndexDeriv[isr][i] * nz_dat + _zIndex[isr][j]] += _xWeightDeriv[isr][i] * _zWeight[isr][j] * pmod[isr * nz_mod + it];
pdat[_xIndex[isr][j] * nz_dat + _zIndexDeriv[isr][i]] += _xWeight[isr][j] * _zWeightDeriv[isr][i] * pmod[(_nsr+isr) * nz_mod + it];
}
}
}
}
......@@ -533,7 +639,7 @@ void injecter_delta_derivative_m5::findIndexDeriv(const std::shared_ptr<hypercub
if (_xIndexDeriv[isr][2] <= 1) _xIndexDeriv[isr][2] = 2;
if (_xIndexDeriv[isr][1] >= nx - 2) _xIndexDeriv[isr][1] = nx - 3;
if (_zIndexDeriv[isr][2] <= 1) _zIndexDeriv[isr][2] = 2;
if (_zIndexDeriv[isr][1] >= nx - 2) _zIndexDeriv[isr][1] = nx - 3;
if (_zIndexDeriv[isr][1] >= nz - 2) _zIndexDeriv[isr][1] = nz - 3;
_xIndexDeriv[isr][0] = _xIndexDeriv[isr][2] - 2;
_xIndexDeriv[isr][1] = _xIndexDeriv[isr][2] - 1;
......@@ -594,8 +700,98 @@ void injecter_delta_derivative_m5::forward(bool add, const std::shared_ptr<vec2D
for (int isr = 0; isr < _nsr; isr++){
for (int i=0; i<5; i++){
for (int j=0; j<4; j++){
for (int j=0; j<4; j++){
for (int i=0; i<5; i++){
pdat[_xIndexDeriv[isr][i] * nz_dat + _zIndex[isr][j]] += _xWeightDeriv[isr][i] * _zWeight[isr][j] * pmod[isr * nz_mod + it];
pdat[_xIndex[isr][j] * nz_dat + _zIndexDeriv[isr][i]] += _xWeight[isr][j] * _zWeightDeriv[isr][i] * pmod[(_nsr+isr) * nz_mod + it];
}
}
}
}
void injecter_delta_derivative_m3_s1::findIndexDeriv(const std::shared_ptr<hypercube> range){
data_t dx = range->getAxis(2).d;
data_t dz = range->getAxis(1).d;
int nx = range->getAxis(2).n;
int nz = range->getAxis(1).n;
int ox = range->getAxis(2).o;
int oz = range->getAxis(1).o;
for (int isr = 0; isr < _nsr; isr++){
assert( (_xloc[isr] >= ox) && (_xloc[isr] <= ox+(nx-1)*dx) && (_zloc[isr] >= oz) && (_zloc[isr] <= oz+(nz-1)*dz) );
// compute the index of the pivot xi
_xIndexDeriv[isr][1] = std::floor((_xloc[isr] - ox)/dx);
_zIndexDeriv[isr][1] = std::floor((_zloc[isr] - oz)/dz);
// check if the pivot falls at or near the boundary
// if so, shift it forward or backward one or two samples to make room for the other points xi-1 xi+1 and xi+2
if (_xIndexDeriv[isr][1] == 0) _xIndexDeriv[isr][1] = 1;
if (_xIndexDeriv[isr][1] >= nx - 2) _xIndexDeriv[isr][1] = nx - 3;
if (_zIndexDeriv[isr][1] == 0) _zIndexDeriv[isr][1] = 1;
if (_zIndexDeriv[isr][1] >= nz - 2) _zIndexDeriv[isr][1] = nz - 3;
_xIndexDeriv[isr][0] = _xIndexDeriv[isr][1] - 1;
_xIndexDeriv[isr][2] = _xIndexDeriv[isr][1] + 1;
_xIndexDeriv[isr][3] = _xIndexDeriv[isr][1] + 2;
_zIndexDeriv[isr][0] = _zIndexDeriv[isr][1] - 1;
_zIndexDeriv[isr][2] = _zIndexDeriv[isr][1] + 1;
_zIndexDeriv[isr][3] = _zIndexDeriv[isr][1] + 2;
}
}
void injecter_delta_derivative_m3_s1::findWeightDeriv(const HxHz * H){
data_t dx = _range->getAxis(2).d;
data_t dz = _range->getAxis(1).d;
int ox = _range->getAxis(2).o;
int oz = _range->getAxis(1).o;
data_t alpha_x, alpha_z;
for (int isr = 0; isr < _nsr; isr++){
// compute alpha based on the pivot: x = xi + alpha * h ; xi is the pivot
alpha_x = (_xloc[isr] - dx * _xIndexDeriv[isr][1] - ox) / dx;
alpha_z = (_zloc[isr] - dz * _zIndexDeriv[isr][1] - oz) / dz;
// compute the weights based on the formulas given in the hpp file
_xWeightDeriv[isr][0] = (1-alpha_x) / (H->getWeightx(_xIndexDeriv[isr][0]) * dx);
_xWeightDeriv[isr][1] = (-1+2*alpha_x) / (2*H->getWeightx(_xIndexDeriv[isr][1]) * dx);
_xWeightDeriv[isr][2] = (-1+alpha_x) / (H->getWeightx(_xIndexDeriv[isr][2]) * dx);
_xWeightDeriv[isr][3] = (1-2*alpha_x) / (2*H->getWeightx(_xIndexDeriv[isr][3]) * dx);
_zWeightDeriv[isr][0] = (1-alpha_z) / (H->getWeightz(_zIndexDeriv[isr][0]) * dz);
_zWeightDeriv[isr][1] = (-1+2*alpha_z) / (2*H->getWeightz(_zIndexDeriv[isr][1]) * dz);
_zWeightDeriv[isr][2] = (-1+alpha_z) / (H->getWeightz(_zIndexDeriv[isr][2]) * dz);
_zWeightDeriv[isr][3] = (1-2*alpha_z) / (2*H->getWeightz(_zIndexDeriv[isr][3]) * dz);
}
}
void injecter_delta_derivative_m3_s1::forward(bool add, const std::shared_ptr<vec2DReg> mod, std::shared_ptr<vec2DReg> dat) const{
if (add == false) dat->zero();
int it = this->getIt();
int nz_mod = mod->getHyper()->getAxis(1).n;
int nz_dat = dat->getHyper()->getAxis(1).n;
const data_t * pmod = mod->getCVals();
data_t * pdat = dat->getVals();
for (int isr = 0; isr < _nsr; isr++){
for (int j=0; j<2; j++){
for (int i=0; i<4; i++){
pdat[_xIndexDeriv[isr][i] * nz_dat + _zIndex[isr][j]] += _xWeightDeriv[isr][i] * _zWeight[isr][j] * pmod[isr * nz_mod + it];
pdat[_xIndex[isr][j] * nz_dat + _zIndexDeriv[isr][i]] += _xWeight[isr][j] * _zWeightDeriv[isr][i] * pmod[(_nsr+isr) * nz_mod + it];
}
......
src/include/scalar/initializer.hpp
+ 5
- 0
View file @ 66a48942
......@@ -331,6 +331,11 @@ public:
inj = new injecter_delta_m2(xloc, zloc, src->getHyper(), vel->getHyper(), H);
}
else if(type == "delta_m3"){
assert(H != nullptr);
inj = new injecter_delta_m3(xloc, zloc, src->getHyper(), vel->getHyper(), H);
}
else if(type == "delta_m4"){
assert(H != nullptr);
inj = new injecter_delta_m4(xloc, zloc, src->getHyper(), vel->getHyper(), H);
......
src/include/scalar/injecter.hpp
+ 230
- 15
View file @ 66a48942
......@@ -262,6 +262,86 @@ public:
};
// injecter/extractor by approximating the Dirac distribution using H norm operator and 3 moment conditions
// in 1D, 3 points are required to approximate delta:
// di-1 = (alpha^2 - alpha)/(2*h*wi-1)
// di = (1 - alpha^2)/(h*wi)
// di+1 = (alpha^2 + alpha)/(2*h*wi+1)
// alpha=(x-xi)/h where the singularity is @ x ; h is the sampling ; wi-1, wi, wi+1 are the coefficients from H (not including the sampling)
// in 2D, it is the product of the 1D approximations
//
// Near the boundaries, the pivot xi may be shifted one sample to make room for 3 pts ; alpha is computed accordingly
// Note that the adjoint is relative to the H norm, it is not the transpose of the forward
class injecter_delta_m3 : public injecter {
protected:
std::vector<std::vector<int> > _xIndex, _zIndex;
std::vector<std::vector<data_t> > _xWeight, _zWeight; // weights for the forward operator
std::vector<std::vector<data_t> > _xWeight2, _zWeight2; // weights for the adjoint operator
public:
// default constructor
injecter_delta_m3(){}
// destructor
~injecter_delta_m3(){};
// constructor
injecter_delta_m3(const std::vector<data_t> &xloc, const std::vector<data_t> &zloc,
const std::shared_ptr<hypercube> domain, const std::shared_ptr<hypercube> range) : injecter(xloc,zloc,domain,range){
_xIndex.resize(_nsr, std::vector<int>(3,0));
_zIndex.resize(_nsr, std::vector<int>(3,0));
_xWeight.resize(_nsr, std::vector<data_t>(3,0.0));
_zWeight.resize(_nsr, std::vector<data_t>(3,0.0));
_xWeight2.resize(_nsr, std::vector<data_t>(3,0.0));
_zWeight2.resize(_nsr, std::vector<data_t>(3,0.0));
this->findIndex(range);
}
// another constructor with the H norm
injecter_delta_m3(const std::vector<data_t> &xloc, const std::vector<data_t> &zloc,
const std::shared_ptr<hypercube> domain, const std::shared_ptr<hypercube> range, const HxHz * H) : injecter(xloc,zloc,domain,range){
_xIndex.resize(_nsr, std::vector<int>(3,0));
_zIndex.resize(_nsr, std::vector<int>(3,0));
_xWeight.resize(_nsr, std::vector<data_t>(3,0.0));
_zWeight.resize(_nsr, std::vector<data_t>(3,0.0));
_xWeight2.resize(_nsr, std::vector<data_t>(3,0.0));
_zWeight2.resize(_nsr, std::vector<data_t>(3,0.0));
this->findIndex(range);
H->getDomain()->checkSame(_range);
this->findWeight(H);
}
// clone the object
injecter_delta_m3 * clone() const {
injecter_delta_m3 * inj = new injecter_delta_m3(_xloc, _zloc, _domain, _range);
inj->_xIndex = _xIndex;
inj->_zIndex = _zIndex;
inj->_xWeight = _xWeight;
inj->_zWeight = _zWeight;
inj->_xWeight2 = _xWeight2;
inj->_zWeight2 = _zWeight2;
return inj;
}
// find the nearest 3 indices of sources/receivers in a given wavefield in each direction
void findIndex(const std::shared_ptr<hypercube> range);
// find the appropriate weights from the H norm
void findWeight(const HxHz * H);
virtual void forward(bool add, const std::shared_ptr<vec2DReg> mod, std::shared_ptr<vec2DReg> dat) const;
virtual void adjoint(bool add, std::shared_ptr<vec2DReg> mod, const std::shared_ptr<vec2DReg> dat) const;
};
// injecter/extractor by approximating the Dirac distribution using H norm operator and 4 moment conditions
// in 1D, 4 points are required to approximate delta:
// di-1 = (-alpha^3 + 3*alpha^2 - 2*alpha)/(6*h*wi-1)
......@@ -354,7 +434,7 @@ public:
// the first time function will be injected with delta derivative in x and normal delta in z
// the second time function will be injected with delta derivative in z and normal delta in x
class injecter_delta_derivative_m3 : public injecter_delta_m2 {
class injecter_delta_derivative_m3 : public injecter_delta_m3 {
protected:
std::vector<std::vector<int> > _xIndexDeriv, _zIndexDeriv;
......@@ -385,16 +465,22 @@ public:
_domain = domain->clone();
_range = range->clone();
_xIndex.resize(_nsr, std::vector<int>(2,0));
_zIndex.resize(_nsr, std::vector<int>(2,0));
//_xIndex.resize(_nsr, std::vector<int>(2,0));
//_zIndex.resize(_nsr, std::vector<int>(2,0));
_xIndex.resize(_nsr, std::vector<int>(3,0));
_zIndex.resize(_nsr, std::vector<int>(3,0));
_xIndexDeriv.resize(_nsr, std::vector<int>(3,0));
_zIndexDeriv.resize(_nsr, std::vector<int>(3,0));
_xWeight.resize(_nsr, std::vector<data_t>(2,0.0));
_zWeight.resize(_nsr, std::vector<data_t>(2,0.0));
//_xWeight.resize(_nsr, std::vector<data_t>(2,0.0));
//_zWeight.resize(_nsr, std::vector<data_t>(2,0.0));
_xWeight.resize(_nsr, std::vector<data_t>(3,0.0));
_zWeight.resize(_nsr, std::vector<data_t>(3,0.0));
_xWeightDeriv.resize(_nsr, std::vector<data_t>(3,0.0));
_zWeightDeriv.resize(_nsr, std::vector<data_t>(3,0.0));
_xWeight2.resize(_nsr, std::vector<data_t>(2,0.0));
_zWeight2.resize(_nsr, std::vector<data_t>(2,0.0));
//_xWeight2.resize(_nsr, std::vector<data_t>(2,0.0));
//_zWeight2.resize(_nsr, std::vector<data_t>(2,0.0));
_xWeight2.resize(_nsr, std::vector<data_t>(3,0.0));
_zWeight2.resize(_nsr, std::vector<data_t>(3,0.0));
this->findIndex(range);
this->findIndexDeriv(range);
......@@ -418,22 +504,29 @@ public:
_domain = domain->clone();
_range = range->clone();
_xIndex.resize(_nsr, std::vector<int>(2,0));
_zIndex.resize(_nsr, std::vector<int>(2,0));
//_xIndex.resize(_nsr, std::vector<int>(2,0));
//_zIndex.resize(_nsr, std::vector<int>(2,0));
_xIndex.resize(_nsr, std::vector<int>(3,0));
_zIndex.resize(_nsr, std::vector<int>(3,0));
_xIndexDeriv.resize(_nsr, std::vector<int>(3,0));
_zIndexDeriv.resize(_nsr, std::vector<int>(3,0));
_xWeight.resize(_nsr, std::vector<data_t>(2,0.0));
_zWeight.resize(_nsr, std::vector<data_t>(2,0.0));
//_xWeight.resize(_nsr, std::vector<data_t>(2,0.0));
//_zWeight.resize(_nsr, std::vector<data_t>(2,0.0));
_xWeight.resize(_nsr, std::vector<data_t>(3,0.0));
_zWeight.resize(_nsr, std::vector<data_t>(3,0.0));
_xWeightDeriv.resize(_nsr, std::vector<data_t>(3,0.0));
_zWeightDeriv.resize(_nsr, std::vector<data_t>(3,0.0));
_xWeight2.resize(_nsr, std::vector<data_t>(2,0.0));
_zWeight2.resize(_nsr, std::vector<data_t>(2,0.0));
//_xWeight2.resize(_nsr, std::vector<data_t>(2,0.0));
//_zWeight2.resize(_nsr, std::vector<data_t>(2,0.0));
_xWeight2.resize(_nsr, std::vector<data_t>(3,0.0));
_zWeight2.resize(_nsr, std::vector<data_t>(3,0.0));
this->findIndex(range);
this->findIndexDeriv(range);
H->getDomain()->checkSame(_range);
this->findWeight(H);
this->findWeightDeriv(H);
}
// clone the object
......@@ -461,9 +554,11 @@ public:
void findWeightDeriv(const HxHz * H);
void forward(bool add, const std::shared_ptr<vec2DReg> mod, std::shared_ptr<vec2DReg> dat) const;
//void adjoint(bool add, std::shared_ptr<vec2DReg> mod, const std::shared_ptr<vec2DReg> dat) const;
};
// injecter/extractor by approximating the Dirac distribution and its spatial derivative using H norm operator and 4-5 moment conditions
// this class is used to inject moment tensor ; for each source, 2 time functions are required
// in 1D, 5 points are required to approximate the derivative of delta:
......@@ -553,6 +648,7 @@ public:
H->getDomain()->checkSame(_range);
this->findWeight(H);
this->findWeightDeriv(H);
}
// clone the object
......@@ -580,6 +676,125 @@ public:
void findWeightDeriv(const HxHz * H);
void forward(bool add, const std::shared_ptr<vec2DReg> mod, std::shared_ptr<vec2DReg> dat) const;
//void adjoint(bool add, std::shared_ptr<vec2DReg> mod, const std::shared_ptr<vec2DReg> dat) const;
};
// injecter/extractor by approximating the Dirac distribution and its spatial derivative using H norm operator and 2-3 moment conditions
// one smoothness condition for the delta derivative
// this class is used to inject moment tensor ; for each source, 2 time functions are required
// in 1D, 4 points are required to approximate the derivative of delta:
// di-1= (1-alpha)/(h*wi-1) ; di=(2*alpha-1)/(2*h*wi) ; di+1=(alpha-1)/(h*wi+1) ; di+2=(-2*alpha+1)/(2*h*wi+2)
// alpha=(x-xi)/h where the singularity is @ x ; h is the sampling ; wi-1, wi, wi+1 & wi+2 are the coefficients from H (including the sampling)
// the first time function will be injected with delta derivative in x and normal delta in z
// the second time function will be injected with delta derivative in z and normal delta in x
class injecter_delta_derivative_m3_s1 : public injecter_delta_m3 {
protected:
std::vector<std::vector<int> > _xIndexDeriv, _zIndexDeriv;
std::vector<std::vector<data_t> > _xWeightDeriv, _zWeightDeriv; // weights for the forward operator
public:
// default constructor
injecter_delta_derivative_m3_s1(){}
// destructor
~injecter_delta_derivative_m3_s1(){};
// constructor
injecter_delta_derivative_m3_s1(const std::vector<data_t> &xloc, const std::vector<data_t> &zloc,
const std::shared_ptr<hypercube> domain, const std::shared_ptr<hypercube> range){
if (xloc.size() != zloc.size())
throw std::logic_error("The xz location vectors doesn't have the same size.\n");
if (domain->getAxis(2).n != 2*xloc.size())
throw std::logic_error("The domain has too many or too little traces.\n");
_it = 0;
_xloc = xloc;
_zloc = zloc;
_nsr = xloc.size();
_domain = domain->clone();
_range = range->clone();
_xIndex.resize(_nsr, std::vector<int>(3,0));
_zIndex.resize(_nsr, std::vector<int>(3,0));
_xIndexDeriv.resize(_nsr, std::vector<int>(4,0));
_zIndexDeriv.resize(_nsr, std::vector<int>(4,0));
_xWeight.resize(_nsr, std::vector<data_t>(3,0.0));
_zWeight.resize(_nsr, std::vector<data_t>(3,0.0));
_xWeightDeriv.resize(_nsr, std::vector<data_t>(4,0.0));
_zWeightDeriv.resize(_nsr, std::vector<data_t>(4,0.0));
_xWeight2.resize(_nsr, std::vector<data_t>(3,0.0));
_zWeight2.resize(_nsr, std::vector<data_t>(3,0.0));
this->findIndex(range);
this->findIndexDeriv(range);
}
// another constructor with the H norm
injecter_delta_derivative_m3_s1(const std::vector<data_t> &xloc, const std::vector<data_t> &zloc,
const std::shared_ptr<hypercube> domain, const std::shared_ptr<hypercube> range, const HxHz * H){
if (xloc.size() != zloc.size())
throw std::logic_error("The xz location vectors doesn't have the same size.\n");
if (domain->getAxis(2).n != 2*xloc.size())
throw std::logic_error("The domain has too many or too little traces.\n");
_it = 0;
_xloc = xloc;
_zloc = zloc;
_nsr = xloc.size();
_domain = domain->clone();
_range = range->clone();
_xIndex.resize(_nsr, std::vector<int>(3,0));
_zIndex.resize(_nsr, std::vector<int>(3,0));
_xIndexDeriv.resize(_nsr, std::vector<int>(4,0));
_zIndexDeriv.resize(_nsr, std::vector<int>(4,0));
_xWeight.resize(_nsr, std::vector<data_t>(3,0.0));
_zWeight.resize(_nsr, std::vector<data_t>(3,0.0));
_xWeightDeriv.resize(_nsr, std::vector<data_t>(4,0.0));
_zWeightDeriv.resize(_nsr, std::vector<data_t>(4,0.0));
_xWeight2.resize(_nsr, std::vector<data_t>(3,0.0));
_zWeight2.resize(_nsr, std::vector<data_t>(3,0.0));
this->findIndex(range);
this->findIndexDeriv(range);
H->getDomain()->checkSame(_range);
this->findWeight(H);
this->findWeightDeriv(H);
}
// clone the object
injecter_delta_derivative_m3_s1 * clone() const {
injecter_delta_derivative_m3_s1 * inj = new injecter_delta_derivative_m3_s1(_xloc, _zloc, _domain, _range);
inj->_xIndex = _xIndex;
inj->_zIndex = _zIndex;
inj->_xIndexDeriv = _xIndexDeriv;
inj->_zIndexDeriv = _zIndexDeriv;
inj->_xWeight = _xWeight;
inj->_zWeight = _zWeight;
inj->_xWeightDeriv = _xWeightDeriv;
inj->_zWeightDeriv = _zWeightDeriv;
inj->_xWeight2 = _xWeight2;
inj->_zWeight2 = _zWeight2;
return inj;
}
// find the nearest 4 indices of sources/receivers in a given wavefield in each direction
void findIndexDeriv(const std::shared_ptr<hypercube> range);
// find the appropriate weights from the H norm
void findWeightDeriv(const HxHz * H);
void forward(bool add, const std::shared_ptr<vec2DReg> mod, std::shared_ptr<vec2DReg> dat) const;
};
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment