zhangyiss/femtic
1
0
Fork 0
mirror of https://github.com/yoshiya-usui/femtic.git synced 2025-05-12 19:51:12 +08:00
Code Issues Actions Packages Projects Releases Wiki Activity
v4.1
femtic/src/ObservedDataStationNMT.cpp

613 lines
26 KiB
C++
Raw Permalink Normal View History

add new file
2021-11-09 00:06:52 +08:00
//-------------------------------------------------------------------------------------------------------
// The MIT License (MIT)
//
// Copyright (c) 2021 Yoshiya Usui
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
//-------------------------------------------------------------------------------------------------------
#include <fstream>
#include <iostream>
#include <string>
#include <algorithm>
#include <assert.h>
#include <iomanip>
#include "ObservedDataStationNMT.h"
#include "OutputFiles.h"
#include "CommonParameters.h"
#include "ResistivityBlock.h"
// Constructer
ObservedDataStationNMT::ObservedDataStationNMT():
ObservedDataStation(),
m_YxObserved(NULL),
m_YyObserved(NULL),
m_YxSD(NULL),
m_YySD(NULL),
m_YxCalculated(NULL),
m_YyCalculated(NULL),
m_YxResidual(NULL),
m_YyResidual(NULL),
m_numElementsIncludingDipole(0),
m_elementsIncludingDipole(NULL),
m_facesIncludingDipole(NULL),
m_localCoordinateValuesStartPoint(NULL),
m_localCoordinateValuesEndPoint(NULL),
m_areaCoordinateValuesStartPoint(NULL),
m_areaCoordinateValuesEndPoint(NULL),
m_rhsVectorIDOfVoltageDifference(0),
m_dataIDOfYx(NULL),
m_dataIDOfYy(NULL)
{
m_location.startPoint.X = 0.0;
m_location.startPoint.Y = 0.0;
m_location.endPoint.X = 0.0;
m_location.endPoint.Y = 0.0;
for( int i = 0; i < 2; ++i ){
m_voltageCalculated[i] = std::complex<double>(0.0,0.0);
}
}
// Destructer
ObservedDataStationNMT::~ObservedDataStationNMT(){
if( m_YxObserved != NULL){
delete[] m_YxObserved;
m_YxObserved = NULL;
}
if( m_YyObserved != NULL){
delete[] m_YyObserved;
m_YyObserved = NULL;
}
if( m_YxSD != NULL){
delete[] m_YxSD;
m_YxSD = NULL;
}
if( m_YySD != NULL){
delete[] m_YySD;
m_YySD = NULL;
}
if( m_YxCalculated != NULL){
delete[] m_YxCalculated;
m_YxCalculated = NULL;
}
if( m_YyCalculated != NULL){
delete[] m_YyCalculated;
m_YyCalculated = NULL;
}
if( m_YxResidual != NULL){
delete[] m_YxResidual;
m_YxResidual = NULL;
}
if( m_YyResidual != NULL){
delete[] m_YyResidual;
m_YyResidual = NULL;
}
if( m_elementsIncludingDipole != NULL){
delete[] m_elementsIncludingDipole;
m_elementsIncludingDipole = NULL;
}
if( m_facesIncludingDipole != NULL){
delete[] m_facesIncludingDipole;
m_facesIncludingDipole = NULL;
}
if( m_localCoordinateValuesStartPoint != NULL){
delete[] m_localCoordinateValuesStartPoint;
m_localCoordinateValuesStartPoint = NULL;
}
if( m_localCoordinateValuesEndPoint != NULL){
delete[] m_localCoordinateValuesEndPoint;
m_localCoordinateValuesEndPoint = NULL;
}
if( m_areaCoordinateValuesStartPoint != NULL){
delete[] m_areaCoordinateValuesStartPoint;
m_areaCoordinateValuesStartPoint = NULL;
}
if( m_areaCoordinateValuesEndPoint != NULL){
delete[] m_areaCoordinateValuesEndPoint;
m_areaCoordinateValuesEndPoint = NULL;
}
if( m_dataIDOfYx != NULL){
delete[] m_dataIDOfYx;
m_dataIDOfYx = NULL;
}
if( m_dataIDOfYy != NULL){
delete[] m_dataIDOfYy;
m_dataIDOfYy = NULL;
}
}
// Read data from input file
void ObservedDataStationNMT::inputObservedData( std::ifstream& inFile ){
inFile >> m_stationID;
inFile >> m_IDOfMagneticFieldStation;
OutputFiles::m_logFile << "# " << std::setw(15) << std::left << m_stationID << std::setw(15) << std::left << m_IDOfMagneticFieldStation << std::endl;
double dbuf(0.0);
inFile >> dbuf;
m_location.startPoint.X = dbuf * CommonParameters::convKilometerToMeter;
inFile >> dbuf;
m_location.startPoint.Y = dbuf * CommonParameters::convKilometerToMeter;
inFile >> dbuf;
m_location.endPoint.X = dbuf * CommonParameters::convKilometerToMeter;
inFile >> dbuf;
m_location.endPoint.Y = dbuf * CommonParameters::convKilometerToMeter;
inFile >> m_numOfFrequency;
const int nFreq = m_numOfFrequency;
if( nFreq > 0 ){
m_freq = new double[nFreq];
m_YxObserved = new std::complex<double>[nFreq];
m_YyObserved = new std::complex<double>[nFreq];
m_YxSD = new CommonParameters::DoubleComplexValues[nFreq];
m_YySD = new CommonParameters::DoubleComplexValues[nFreq];
//m_YxCalculated = new std::complex<double>[nFreq];
//m_YyCalculated = new std::complex<double>[nFreq];
for( int i = 0; i < nFreq; ++i ){
inFile >> m_freq[i];
double dbuf1(0.0);
double dbuf2(0.0);
inFile >> dbuf1 >> dbuf2;
m_YxObserved[i] = std::complex<double>(dbuf1,dbuf2);
inFile >> dbuf1 >> dbuf2;
m_YyObserved[i] = std::complex<double>(dbuf1,dbuf2);
inFile >> m_YxSD[i].realPart;
inFile >> m_YxSD[i].imagPart;
inFile >> m_YySD[i].realPart;
inFile >> m_YySD[i].imagPart;
}
}
#ifdef _DEBUG_WRITE
std::cout << " NMT " << m_stationID << " " << m_IDOfMagneticFieldStation << std::endl;
std::cout << m_location.startPoint.X << " "
<< m_location.startPoint.Y << " "
<< m_location.endPoint.X << " "
<< m_location.endPoint.Y << std::endl;
std::cout << m_numOfFrequency << std::endl;
for( int i = 0; i < m_numOfFrequency; ++i ){
std::cout << m_freq[i] << " "
<< m_YxObserved[i] << " "
<< m_YyObserved[i] << " "
<< m_YxSD[i].realPart << " "
<< m_YxSD[i].imagPart << " "
<< m_YySD[i].realPart << " "
<< m_YySD[i].imagPart << std::endl;
}
#endif
}
// Find elements including dipole
void ObservedDataStationNMT::findElementsIncludingDipole(){
if( ( AnalysisControl::getInstance() )->getTypeOfMesh() == MeshData::TETRA ){// Tetra/Prism mesh
const MeshDataTetraElement* const ptrMeshDataTetraElement = ( AnalysisControl::getInstance() )->getPointerOfMeshDataTetraElement();
std::vector<CommonParameters::AreaCoords> localCoordStartPoint;
std::vector<CommonParameters::AreaCoords> localCoordEndPoint;
std::vector<int> elementsIncludingDipole;
std::vector<int> facesIncludingDipole;
ptrMeshDataTetraElement->findElementsIncludingDipoleOnSurface( m_location.startPoint.X, m_location.startPoint.Y, m_location.endPoint.X, m_location.endPoint.Y,
elementsIncludingDipole, facesIncludingDipole, localCoordStartPoint, localCoordEndPoint );
m_numElementsIncludingDipole = static_cast<int>( elementsIncludingDipole.size() );
if( m_numElementsIncludingDipole > 0 ){
m_elementsIncludingDipole = new int[m_numElementsIncludingDipole];
m_facesIncludingDipole = new int[m_numElementsIncludingDipole];
m_areaCoordinateValuesStartPoint = new CommonParameters::AreaCoords[m_numElementsIncludingDipole];
m_areaCoordinateValuesEndPoint = new CommonParameters::AreaCoords[m_numElementsIncludingDipole];
}
for( int i = 0; i < m_numElementsIncludingDipole; ++i ){
m_elementsIncludingDipole[i] = elementsIncludingDipole[i];
m_facesIncludingDipole[i] = facesIncludingDipole[i];
m_areaCoordinateValuesStartPoint[i] = localCoordStartPoint[i];
m_areaCoordinateValuesEndPoint[i] = localCoordEndPoint[i];
}
#ifdef _DEBUG_WRITE
std::cout << "m_stationID " << m_stationID << std::endl;
for( int i = 0; i < m_numElementsIncludingDipole; ++i ){
std::cout << m_elementsIncludingDipole[i] << " " << m_facesIncludingDipole[i] << " " << m_areaCoordinateValuesStartPoint[i].coord0 << " " << m_areaCoordinateValuesStartPoint[i].coord1 << " " << m_areaCoordinateValuesStartPoint[i].coord2 << std::endl;
std::cout << m_elementsIncludingDipole[i] << " " << m_facesIncludingDipole[i] << " " << m_areaCoordinateValuesStartPoint[i].coord0 << " " << m_areaCoordinateValuesStartPoint[i].coord1 << " " << m_areaCoordinateValuesStartPoint[i].coord2 << std::endl;
}
#endif
}else{// Hexa mesh
const MeshDataBrickElement* const ptrMeshDataBrickElement = ( AnalysisControl::getInstance() )->getPointerOfMeshDataBrickElement();
std::vector<double> localCoordXStartPoint;
std::vector<double> localCoordYStartPoint;
std::vector<double> localCoordXEndPoint;
std::vector<double> localCoordYEndPoint;
std::vector<int> elementsIncludingDipole;
ptrMeshDataBrickElement->findElementsIncludingDipoleOnSurface( m_location.startPoint.X, m_location.startPoint.Y, m_location.endPoint.X, m_location.endPoint.Y,
elementsIncludingDipole, localCoordXStartPoint, localCoordYStartPoint, localCoordXEndPoint, localCoordYEndPoint );
m_numElementsIncludingDipole = static_cast<int>( elementsIncludingDipole.size() );
if( m_numElementsIncludingDipole > 0 ){
m_elementsIncludingDipole = new int[m_numElementsIncludingDipole];
m_localCoordinateValuesStartPoint = new CommonParameters::locationXY[m_numElementsIncludingDipole];
m_localCoordinateValuesEndPoint = new CommonParameters::locationXY[m_numElementsIncludingDipole];
}
for( int i = 0; i < m_numElementsIncludingDipole; ++i ){
m_elementsIncludingDipole[i] = elementsIncludingDipole[i];
m_localCoordinateValuesStartPoint[i].X = localCoordXStartPoint[i];
m_localCoordinateValuesStartPoint[i].Y = localCoordYStartPoint[i];
m_localCoordinateValuesEndPoint[i].X = localCoordXEndPoint[i];
m_localCoordinateValuesEndPoint[i].Y = localCoordYEndPoint[i];
}
#ifdef _DEBUG_WRITE
std::cout << "m_stationID " << m_stationID << std::endl;
std::cout << "m_elementsIncludingDipole m_localCoordinateValuesStartPoint[i].x m_localCoordinateValuesStartPoint[i].y m_localCoordinateValuesEndPoint[i].x m_localCoordinateValuesEndPoint[i].y" << std::endl;
for( int i = 0; i < m_numElementsIncludingDipole; ++i ){
std::cout << m_elementsIncludingDipole[i] << " " << m_localCoordinateValuesStartPoint[i].X << " " << m_localCoordinateValuesStartPoint[i].Y << " " << m_localCoordinateValuesEndPoint[i].X << " " << m_localCoordinateValuesEndPoint[i].Y << std::endl;
}
#endif
}
}
// Calulate difference of voltages
void ObservedDataStationNMT::calculateVoltageDifferences( const Forward3D* const ptrForward3D, const int rhsVectorIDOfVoltageDifference ){
const int iPol = ptrForward3D->getPolarizationCurrent();
if( ( AnalysisControl::getInstance() )->getTypeOfMesh() == MeshData::TETRA ){// Tetra mesh
m_voltageCalculated[iPol] = ptrForward3D->calcVoltageDifference( m_numElementsIncludingDipole, m_elementsIncludingDipole, m_facesIncludingDipole, m_areaCoordinateValuesStartPoint, m_areaCoordinateValuesEndPoint );
}else{// Hexa mesh
m_voltageCalculated[iPol] = ptrForward3D->calcVoltageDifference( m_numElementsIncludingDipole, m_elementsIncludingDipole, m_localCoordinateValuesStartPoint, m_localCoordinateValuesEndPoint );
}
#ifdef _DEBUG_WRITE
std::cout << "iPol V : " << iPol << " " << m_voltageCalculated[iPol] << std::endl;
#endif
// For inversion
m_rhsVectorIDOfVoltageDifference = rhsVectorIDOfVoltageDifference;
}
// Calulate Network-MT response
//void ObservedDataStationNMT::calculateNetworkMTResponse( const int freqIDAmongThisPE, const ObservedDataStationPoint* const ptrStationOfMagneticField, int& icount ){
void ObservedDataStationNMT::calculateNetworkMTResponse( const double freq, const ObservedDataStationPoint* const ptrStationOfMagneticField, int& icount ){
const int freqIDThisPEInSta = getFreqIDsAmongThisPE( freq );
if( freqIDThisPEInSta < 0 ){// Specified frequency is not the ones calculated by this PE in this station
return;
}
const int freqIDGlobalInSta = m_freqIDsAmongThisStationCalculatedByThisPE[ freqIDThisPEInSta ];
const std::complex<double> HxCalculated[2] = { ptrStationOfMagneticField->getHxCalculated( 0 ), ptrStationOfMagneticField->getHxCalculated( 1 ) };
const std::complex<double> HyCalculated[2] = { ptrStationOfMagneticField->getHyCalculated( 0 ), ptrStationOfMagneticField->getHyCalculated( 1 ) };
const std::complex<double> det = HxCalculated[0]*HyCalculated[1] - HxCalculated[1]*HyCalculated[0];
m_YxCalculated[freqIDThisPEInSta] = ( m_voltageCalculated[0]*HyCalculated[1] - m_voltageCalculated[1]*HyCalculated[0] ) / det;
m_YyCalculated[freqIDThisPEInSta] = ( m_voltageCalculated[1]*HxCalculated[0] - m_voltageCalculated[0]*HxCalculated[1] ) / det;
m_YxResidual[freqIDThisPEInSta].realPart = ( m_YxObserved[freqIDGlobalInSta].real() - m_YxCalculated[freqIDThisPEInSta].real() ) / m_YxSD[freqIDGlobalInSta].realPart;
m_YxResidual[freqIDThisPEInSta].imagPart = ( m_YxObserved[freqIDGlobalInSta].imag() - m_YxCalculated[freqIDThisPEInSta].imag() ) / m_YxSD[freqIDGlobalInSta].imagPart;
m_YyResidual[freqIDThisPEInSta].realPart = ( m_YyObserved[freqIDGlobalInSta].real() - m_YyCalculated[freqIDThisPEInSta].real() ) / m_YySD[freqIDGlobalInSta].realPart;
m_YyResidual[freqIDThisPEInSta].imagPart = ( m_YyObserved[freqIDGlobalInSta].imag() - m_YyCalculated[freqIDThisPEInSta].imag() ) / m_YySD[freqIDGlobalInSta].imagPart;
#ifdef _DEBUG_WRITE
std::cout << "ifreq Yx Yy: " << freqIDThisPEInSta << " " << m_YxCalculated[freqIDThisPEInSta] << " " << m_YyCalculated[freqIDThisPEInSta]<< std::endl;
#endif
// For inversion
//ObservedData* const ptrObservedData = ObservedData::getInstance();
//m_dataIDOfYx[freqIDThisPEInSta] = ptrObservedData->incrementNumObservedDataThisPE( freqIDThisPEInSta );
//m_dataIDOfYy[freqIDThisPEInSta] = ptrObservedData->incrementNumObservedDataThisPE( freqIDThisPEInSta );
m_dataIDOfYx[freqIDThisPEInSta].realPart = icount++;
m_dataIDOfYx[freqIDThisPEInSta].imagPart = icount++;
m_dataIDOfYy[freqIDThisPEInSta].realPart = icount++;
m_dataIDOfYy[freqIDThisPEInSta].imagPart = icount++;
}
// Initialize difference of voltages
void ObservedDataStationNMT::initializeVoltageDifferences( const int iPol ){
m_voltageCalculated[iPol] = std::complex<double>(0.0,0.0);
}
// Initialize Network-MT responses and errors
void ObservedDataStationNMT::initializeNetworkMTResponsesAndErrors(){
for( int i = 0; i < m_numOfFreqCalculatedByThisStaAndPE; ++i ){
m_YxResidual[i].realPart = 0.0;
m_YxResidual[i].imagPart = 0.0;
m_YyResidual[i].realPart = 0.0;
m_YyResidual[i].imagPart = 0.0;
m_YxCalculated[i] = std::complex<double>(0.0,0.0);
m_YyCalculated[i] = std::complex<double>(0.0,0.0);
}
}
// Allocate memory for the calculated values of Network-MT responses and errors
void ObservedDataStationNMT::allocateMemoryForCalculatedValues(){
if( m_YxCalculated != NULL){
delete[] m_YxCalculated;
m_YxCalculated = NULL;
}
if( m_YyCalculated != NULL){
delete[] m_YyCalculated;
m_YyCalculated = NULL;
}
if( m_YxResidual != NULL){
delete[] m_YxResidual;
m_YxResidual = NULL;
}
if( m_YyResidual != NULL){
delete[] m_YyResidual;
m_YyResidual = NULL;
}
if( m_dataIDOfYx != NULL){
delete[] m_dataIDOfYx;
m_dataIDOfYx = NULL;
}
if( m_dataIDOfYy != NULL){
delete[] m_dataIDOfYy;
m_dataIDOfYy = NULL;
}
if( m_numOfFreqCalculatedByThisStaAndPE > 0 ){
// If total number of frequencies calculated by this PE is not more than zero, do not allocate memory
m_YxCalculated = new std::complex<double>[m_numOfFreqCalculatedByThisStaAndPE];
m_YyCalculated = new std::complex<double>[m_numOfFreqCalculatedByThisStaAndPE];
m_YxResidual = new CommonParameters::DoubleComplexValues[m_numOfFreqCalculatedByThisStaAndPE];
m_YyResidual = new CommonParameters::DoubleComplexValues[m_numOfFreqCalculatedByThisStaAndPE];
m_dataIDOfYx = new CommonParameters::InitComplexValues[m_numOfFreqCalculatedByThisStaAndPE];
m_dataIDOfYy = new CommonParameters::InitComplexValues[m_numOfFreqCalculatedByThisStaAndPE];
// Initialize
for( int i = 0; i < m_numOfFreqCalculatedByThisStaAndPE; ++i ){
m_YxCalculated[i] = std::complex<double>(0.0,0.0);
m_YyCalculated[i] = std::complex<double>(0.0,0.0);;
m_YxResidual[i].realPart = 0.0;
m_YxResidual[i].imagPart = 0.0;
m_YyResidual[i].realPart = 0.0;
m_YyResidual[i].imagPart = 0.0;
m_dataIDOfYx[i].realPart = -1;
m_dataIDOfYx[i].imagPart = -1;
m_dataIDOfYy[i].realPart = -1;
m_dataIDOfYy[i].imagPart = -1;
}
}
}
// Output calculated values of Network-MT responses
void ObservedDataStationNMT::outputCalculatedValues() const{
int icount(0);
for( std::vector<int>::const_iterator itr = m_freqIDsAmongThisStationCalculatedByThisPE.begin(); itr != m_freqIDsAmongThisStationCalculatedByThisPE.end(); ++itr ){
//OutputFiles::m_csvFile << m_stationID << "," << m_freq[*itr] << "," <<
// m_YxCalculated[icount].real() << "," << m_YxCalculated[icount].imag() << "," <<
// m_YyCalculated[icount].real() << "," << m_YyCalculated[icount].imag() << "," <<
// m_YxResidual[icount].realPart << "," << m_YxResidual[icount].imagPart << "," <<
// m_YyResidual[icount].realPart << "," << m_YyResidual[icount].imagPart << "," <<
// m_YxObserved[*itr].real() << "," << m_YxObserved[*itr].imag() << "," <<
// m_YyObserved[*itr].real() << "," << m_YyObserved[*itr].imag() << "," <<
// m_YxSD[*itr].realPart << "," << m_YxSD[*itr].imagPart << "," <<
// m_YySD[*itr].realPart << "," << m_YySD[*itr].imagPart << std::endl;
fprintf( OutputFiles::m_csvFile, "%10d, %15.6e, %15.6e, %15.6e, %15.6e, %15.6e,",
m_stationID, m_freq[*itr], m_YxCalculated[icount].real(), m_YxCalculated[icount].imag(), m_YyCalculated[icount].real(), m_YyCalculated[icount].imag() );
fprintf( OutputFiles::m_csvFile, " %15.6e, %15.6e, %15.6e, %15.6e,",
m_YxResidual[icount].realPart, m_YxResidual[icount].imagPart, m_YyResidual[icount].realPart, m_YyResidual[icount].imagPart );
fprintf( OutputFiles::m_csvFile, " %15.6e, %15.6e, %15.6e, %15.6e,",
m_YxObserved[*itr].real() , m_YxObserved[*itr].imag(), m_YyObserved[*itr].real(), m_YyObserved[*itr].imag() );
fprintf( OutputFiles::m_csvFile, " %15.6e, %15.6e, %15.6e, %15.6e,\n",
m_YxSD[*itr].realPart, m_YxSD[*itr].imagPart, m_YySD[*itr].realPart, m_YySD[*itr].imagPart );
++icount;
}
}
// Calulate interpolator vector of voltage difference
void ObservedDataStationNMT::calcInterpolatorVectorOfVoltageDifference( Forward3D* const ptrForward3D ){
if( ( AnalysisControl::getInstance() )->getTypeOfMesh() == MeshData::TETRA ){// Tetra mesh
ptrForward3D->calcInterpolatorVectorOfVoltageDifference( m_numElementsIncludingDipole, m_elementsIncludingDipole, m_facesIncludingDipole, m_areaCoordinateValuesStartPoint, m_areaCoordinateValuesEndPoint, m_rhsVectorIDOfVoltageDifference );
}else{// Hexa mesh
ptrForward3D->calcInterpolatorVectorOfVoltageDifference( m_numElementsIncludingDipole, m_elementsIncludingDipole, m_localCoordinateValuesStartPoint, m_localCoordinateValuesEndPoint, m_rhsVectorIDOfVoltageDifference );
}
}
// Calulate sensitivity matrix of Impedance tensors
void ObservedDataStationNMT::calculateSensitivityMatrix( const double freq, const int nModel,
const ObservedDataStationPoint* const ptrStationOfMagneticField,
const std::complex<double>* const derivativesOfEMFieldExPol,
const std::complex<double>* const derivativesOfEMFieldEyPol,
double* const sensitivityMatrix ) const{
const int freqIDThisPEInSta = getFreqIDsAmongThisPE( freq );
if( freqIDThisPEInSta < 0 ){// Specified frequency is not the ones calculated by this PE in this station
return;
}
const int freqIDGlobalInSta = m_freqIDsAmongThisStationCalculatedByThisPE[ freqIDThisPEInSta ];
const std::complex<double> HxCalculated[2] = { ptrStationOfMagneticField->getHxCalculated( 0 ), ptrStationOfMagneticField->getHxCalculated( 1 ) };
const std::complex<double> HyCalculated[2] = { ptrStationOfMagneticField->getHyCalculated( 0 ), ptrStationOfMagneticField->getHyCalculated( 1 ) };
const std::complex<double> det = HxCalculated[0]*HyCalculated[1] - HxCalculated[1]*HyCalculated[0];
const std::complex<double> divDet = std::complex<double>(1.0,0.0) / det;
const std::complex<double> divDet2 = divDet * divDet;
const long long rhsVectorIDOfHx = static_cast<long long>(ptrStationOfMagneticField->getRhsVectorIDOfHx());
const long long rhsVectorIDOfHy = static_cast<long long>(ptrStationOfMagneticField->getRhsVectorIDOfHy());
const long long nBlkNotFixed = static_cast<long long>(( ResistivityBlock::getInstance() )->getNumResistivityBlockNotFixed());
for( long long imdl = 0; imdl < nBlkNotFixed; ++imdl ){
const std::complex<double> work1 = derivativesOfEMFieldExPol[ nBlkNotFixed * rhsVectorIDOfHx + imdl ] * HyCalculated[1]
+ derivativesOfEMFieldEyPol[ nBlkNotFixed * rhsVectorIDOfHy + imdl ] * HxCalculated[0]
- derivativesOfEMFieldEyPol[ nBlkNotFixed * rhsVectorIDOfHx + imdl ] * HyCalculated[0]
- derivativesOfEMFieldExPol[ nBlkNotFixed * rhsVectorIDOfHy + imdl ] * HxCalculated[1];
// dYx/dm
const std::complex<double> workX1 = derivativesOfEMFieldExPol[ nBlkNotFixed * static_cast<long long>(m_rhsVectorIDOfVoltageDifference) + imdl ] * HyCalculated[1]
+ derivativesOfEMFieldEyPol[ nBlkNotFixed * rhsVectorIDOfHy + imdl ] * m_voltageCalculated[0]
- derivativesOfEMFieldEyPol[ nBlkNotFixed * static_cast<long long>(m_rhsVectorIDOfVoltageDifference) + imdl ] * HyCalculated[0]
- derivativesOfEMFieldExPol[ nBlkNotFixed * rhsVectorIDOfHy + imdl ] * m_voltageCalculated[1];
const std::complex<double> workX2 = m_voltageCalculated[0]*HyCalculated[1] - m_voltageCalculated[1]*HyCalculated[0];
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfYx[freqIDThisPEInSta].realPart) + imdl ] = ( workX1 * divDet - work1 * workX2 * divDet2 ).real() / m_YxSD[freqIDGlobalInSta].realPart;
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfYx[freqIDThisPEInSta].imagPart) + imdl ] = ( workX1 * divDet - work1 * workX2 * divDet2 ).imag() / m_YxSD[freqIDGlobalInSta].imagPart;
// dYy/dm
const std::complex<double> workY1 = derivativesOfEMFieldEyPol[ nBlkNotFixed * static_cast<long long>(m_rhsVectorIDOfVoltageDifference) + imdl ] * HxCalculated[0]
+ derivativesOfEMFieldExPol[ nBlkNotFixed * rhsVectorIDOfHx + imdl ] * m_voltageCalculated[1]
- derivativesOfEMFieldExPol[ nBlkNotFixed * static_cast<long long>(m_rhsVectorIDOfVoltageDifference) + imdl ] * HxCalculated[1]
- derivativesOfEMFieldEyPol[ nBlkNotFixed * rhsVectorIDOfHx + imdl ] * m_voltageCalculated[0];
const std::complex<double> workY2 = m_voltageCalculated[1]*HxCalculated[0] - m_voltageCalculated[0]*HxCalculated[1];
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfYy[freqIDThisPEInSta].realPart) + imdl ] = ( workY1 * divDet - work1 * workY2 * divDet2 ).real() / m_YySD[freqIDGlobalInSta].realPart;
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfYy[freqIDThisPEInSta].imagPart) + imdl ] = ( workY1 * divDet - work1 * workY2 * divDet2 ).imag() / m_YySD[freqIDGlobalInSta].imagPart;
}
}
// Calculate data vector of this PE
void ObservedDataStationNMT::calculateResidualVectorOfDataThisPE( const double freq, const int offset, double* vector ) const{
const int freqIDThisPEInSta = getFreqIDsAmongThisPE( freq );
if( freqIDThisPEInSta < 0 ){// Specified frequency is not the ones calculated by this PE in this station
return;
}
vector[ offset + m_dataIDOfYx[freqIDThisPEInSta].realPart ] = m_YxResidual[freqIDThisPEInSta].realPart;
vector[ offset + m_dataIDOfYx[freqIDThisPEInSta].imagPart ] = m_YxResidual[freqIDThisPEInSta].imagPart;
vector[ offset + m_dataIDOfYy[freqIDThisPEInSta].realPart ] = m_YyResidual[freqIDThisPEInSta].realPart;
vector[ offset + m_dataIDOfYy[freqIDThisPEInSta].imagPart ] = m_YyResidual[freqIDThisPEInSta].imagPart;
}
// Calulate sum of square of misfit
double ObservedDataStationNMT::calculateErrorSumOfSquaresThisPE() const{
double misfit(0.0);
for( int ifreq = 0; ifreq < m_numOfFreqCalculatedByThisStaAndPE; ++ifreq ){
misfit += pow( m_YxResidual[ifreq].realPart , 2 );
misfit += pow( m_YxResidual[ifreq].imagPart , 2 );
misfit += pow( m_YyResidual[ifreq].realPart , 2 );
misfit += pow( m_YyResidual[ifreq].imagPart , 2 );
}
return misfit;
}
// Get location of the station
const CommonParameters::locationDipole& ObservedDataStationNMT::getLocationOfStation() const{
return m_location;
}
// Get Z coordinate of the point
double ObservedDataStationNMT::getZCoordOfPoint( const int num ) const{
assert( num == 0 || num == 1 );
int iElem(-1);
int iFace(-1);
CommonParameters::AreaCoords areaCoord = { -1.0, -1.0, -1.0 };
if( ( AnalysisControl::getInstance() )->getTypeOfMesh() == MeshData::TETRA ){// Tetra mesh
const MeshDataTetraElement* const ptrMeshDataTetraElement = ( AnalysisControl::getInstance() )->getPointerOfMeshDataTetraElement();
double dummy(0.0);
if( num == 0 ){
iElem = ptrMeshDataTetraElement->findElementIncludingPointOnSurface( m_location.startPoint.X, m_location.startPoint.Y, iFace, areaCoord, false, false, dummy, dummy );
}else{
iElem = ptrMeshDataTetraElement->findElementIncludingPointOnSurface( m_location.endPoint.X, m_location.endPoint.Y, iFace, areaCoord, false, false, dummy, dummy );
}
return ptrMeshDataTetraElement->calcZCoordOfPointOnFace( iElem, iFace, areaCoord );
}
else{
const MeshDataBrickElement* const ptrMeshDataBrickElement = ( AnalysisControl::getInstance() )->getPointerOfMeshDataBrickElement();
double dummy(0.0);
double localCoordZ(0.0);
if( num == 0 ){
iElem = ptrMeshDataBrickElement->findElementIncludingPointOnSurface( m_location.startPoint.X, m_location.startPoint.Y, dummy, dummy, localCoordZ, false, false, dummy, dummy );
}else{
iElem = ptrMeshDataBrickElement->findElementIncludingPointOnSurface( m_location.endPoint.X, m_location.endPoint.Y, dummy, dummy, localCoordZ, false, false, dummy, dummy );
}
return ptrMeshDataBrickElement->calcGlobalCoordZ( iElem, localCoordZ );
}
}
Reference in New Issue Copy Permalink