From 9c41b3e7e9201203ed05e723c8f6f0fbcffc81b2 Mon Sep 17 00:00:00 2001 From: yoshiya-usui Date: Tue, 22 Feb 2022 11:37:11 +0900 Subject: [PATCH] add e-field integration functions for dhexa --- src/CommonParameters.h | 2 +- src/Forward3DBrickElement0thOrder.cpp | 2 +- ...ward3DNonConformingHexaElement0thOrder.cpp | 286 ++++++++++++++++- ...orward3DNonConformingHexaElement0thOrder.h | 17 + src/MeshData.cpp | 192 +++++++++++- src/MeshData.h | 24 ++ src/MeshDataNonConformingHexaElement.cpp | 293 ++++++++++++++++++ src/MeshDataNonConformingHexaElement.h | 12 +- src/MeshDataTetraElement.cpp | 207 +------------ src/MeshDataTetraElement.h | 20 -- src/ObservedData.cpp | 2 + src/ObservedDataStationNMT.cpp | 53 ++++ src/ObservedDataStationNMT2.cpp | 56 +++- 13 files changed, 928 insertions(+), 238 deletions(-) diff --git a/src/CommonParameters.h b/src/CommonParameters.h index 1221915..ba91bcb 100644 --- a/src/CommonParameters.h +++ b/src/CommonParameters.h @@ -162,7 +162,7 @@ static char programName[]="femtic"; // [MajorVersion#].[MinorVersion#].[Revision#] // x.x.xa -> alpha version // x.x.xb -> beta version -static char versionID[]="4.1.9"; +static char versionID[]="4.2"; } diff --git a/src/Forward3DBrickElement0thOrder.cpp b/src/Forward3DBrickElement0thOrder.cpp index fb71c49..67358c3 100644 --- a/src/Forward3DBrickElement0thOrder.cpp +++ b/src/Forward3DBrickElement0thOrder.cpp @@ -647,7 +647,7 @@ std::complex Forward3DBrickElement0thOrder::calcVoltageDifference( const //----- debug >>>>> #ifdef _DEBUG_WRITE std::cout << "areaWithSign calcValueRotatedElectricFieldZDirection : " << areaWithSign << " " << calcValueRotatedElectricFieldZDirection( elememtsIncludingDipole[ielem], localCoordinateValuesStartPoint[ielem].X, localCoordinateValuesStartPoint[ielem].Y, -1.0 ) << std::endl; - std::cout << "voltageDifference : " << voltageDifference << std::endl; + std::cout << "voltageDifference total : " << voltageDifference << std::endl; #endif //----- debug <<<<< diff --git a/src/Forward3DNonConformingHexaElement0thOrder.cpp b/src/Forward3DNonConformingHexaElement0thOrder.cpp index 0043505..af2fa29 100644 --- a/src/Forward3DNonConformingHexaElement0thOrder.cpp +++ b/src/Forward3DNonConformingHexaElement0thOrder.cpp @@ -421,6 +421,34 @@ std::complex Forward3DNonConformingHexaElement0thOrder::calcValueRotated } +// Calculate X component of electric field only from the edges on the Earth's surface +std::complex Forward3DNonConformingHexaElement0thOrder::calcValueRotatedElectricFieldNormal( const int iElem, const double xLocal, const double yLocal ) const{ + + assert( m_solution != NULL ); + assert( m_IDsLocal2Global != NULL ); + + std::complex val(0.0, 0.0); + + Forward3D::Matrix2x2 jacobMat; + const double detJacob = calc2DJacobianMatrixForEarthSurface( iElem, xLocal, yLocal, jacobMat ); + Forward3D::Matrix2x2 invJacobMat; + const double invDet = 1.0 / detJacob; + invJacobMat.mat11 = jacobMat.mat22 * invDet; + invJacobMat.mat12 = - jacobMat.mat12 * invDet; + invJacobMat.mat21 = - jacobMat.mat21 * invDet; + invJacobMat.mat22 = jacobMat.mat11 * invDet; + + const int iFace = 4; + for( int i = 0; i < 4; ++i ){ + const int iEdge = m_MeshDataNonConformingHexaElement.getEdgeIDLocalFromFaceIDLocal( iFace, i ); + const double length = m_MeshDataNonConformingHexaElement.calcEdgeLengthFromElementAndEdge( iElem, iEdge ); + val += m_solution[ m_IDsLocal2Global[iElem][iEdge] ] * std::complex( get2DShapeFuncRotatedForEarthSurface(0.0, 0.0, i, invJacobMat) * length, 0.0 ); + } + + return val; + +} + // Calculate X component of electric field only from the edges on the Earth's surface std::complex Forward3DNonConformingHexaElement0thOrder::calcValueElectricFieldXDirectionFromEdgesOnEarthSurface( const int iElem, const int iFace, const double uCoord, const double vCoord ) const{ OutputFiles::m_logFile << "Error : " << __FUNCTION__ << " is not implemented" << std::endl; @@ -646,6 +674,38 @@ void Forward3DNonConformingHexaElement0thOrder::calcInterpolatorVectorOfRotatedE } +// Calculate interpolator vector of X component of electric field only from the edges on the Earth's surface +void Forward3DNonConformingHexaElement0thOrder::calcInterpolatorVectorOfRotatedElectricFieldNormal( const int iElem, const double xLocal, const double yLocal, + const int irhs, const std::complex& factor ){ + + assert( m_solution != NULL ); + assert( m_IDsLocal2Global != NULL ); + + Forward3D::Matrix2x2 jacobMat; + const double detJacob = calc2DJacobianMatrixForEarthSurface( iElem, xLocal, yLocal, jacobMat ); + Forward3D::Matrix2x2 invJacobMat; + const double invDet = 1.0 / detJacob; + invJacobMat.mat11 = jacobMat.mat22 * invDet; + invJacobMat.mat12 = - jacobMat.mat12 * invDet; + invJacobMat.mat21 = - jacobMat.mat21 * invDet; + invJacobMat.mat22 = jacobMat.mat11 * invDet; + + const int iPol = getPolarizationCurrent(); + const int iFace = 4; + for( int i = 0; i < 4; ++i ){ + const int iEdge = m_MeshDataNonConformingHexaElement.getEdgeIDLocalFromFaceIDLocal( iFace, i ); + const int irow = m_IDsGlobal2AfterDegenerated[iPol][ m_IDsLocal2Global[iElem][iEdge] ]; + if( irow < 0 ){ + continue; + } + const double length = m_MeshDataNonConformingHexaElement.calcEdgeLengthFromElementAndEdge( iElem, iEdge ); + //val += m_solution[ m_IDsLocal2Global[iElem][iEdge] ] * std::complex( get2DShapeFuncRotatedForEarthSurface(0.0, 0.0, i, invJacobMat) * length, 0.0 ); + const std::complex val = std::complex( get2DShapeFuncRotatedForEarthSurface(0.0, 0.0, i, invJacobMat) * length, 0.0 ) * factor; + addValuesToRhsVectorsByConsideringMPC( irow, irhs, val ); + } + +} + // Calculate interpolator vector of X component of electric field only from the edges on the Earth's surface void Forward3DNonConformingHexaElement0thOrder::calcInterpolatorVectorOfElectricFieldXDirectionFromEdgesOnEarthSurface( const int iElem, const int iFace, const double uCoord, const double vCoord, const int irhs, const std::complex& factor ){ OutputFiles::m_logFile << "Error : " << __FUNCTION__ << " is not implemented" << std::endl; @@ -789,9 +849,57 @@ void Forward3DNonConformingHexaElement0thOrder::calcInterpolatorVectorOfMagnetic } // Calculate interpolator vector of difference of voltage -void Forward3DNonConformingHexaElement0thOrder::calcInterpolatorVectorOfVoltageDifference( const int nElem, const int* elememtsIncludingDipole, const CommonParameters::locationXY* localCoordinateValuesStartPoint, const CommonParameters::locationXY* localCoordinateValuesEndPoint, const int irhs ){ - OutputFiles::m_logFile << "Error : " << __FUNCTION__ << " is not implemented" << std::endl; - exit(1); +void Forward3DNonConformingHexaElement0thOrder::calcInterpolatorVectorOfVoltageDifference( const int nElem, const int* elememtsIncludingDipole, + const CommonParameters::locationXY* localCoordinateValuesStartPoint, const CommonParameters::locationXY* localCoordinateValuesEndPoint, const int irhs ){ + + assert( m_IDsLocal2Global != NULL ); + const int iPol = getPolarizationCurrent(); + assert( m_IDsGlobal2AfterDegenerated[iPol] != NULL ); + + for( int ielem = 0; ielem < nElem; ++ielem ){ + const int elemID = elememtsIncludingDipole[ielem]; + const CommonParameters::locationXY startCoord = localCoordinateValuesStartPoint[ielem]; + const CommonParameters::locationXY endCoord = localCoordinateValuesEndPoint[ielem]; + + bool rotationDirectionPlus = true; + CommonParameters::locationXY sharedPoint = {0.0, 0.0}; + const bool integralXCompFirst = doesIntegralXCompFirst( startCoord, endCoord, rotationDirectionPlus, sharedPoint ); + + const double lengthX = m_MeshDataNonConformingHexaElement.getEdgeLengthX(elemID); + const double lengthY = m_MeshDataNonConformingHexaElement.getEdgeLengthY(elemID); + Forward3D::Matrix3x3 JacobMat; + const double detJacob = calcJacobianMatrix( elemID, 0.0, 0.0, -1.0, JacobMat ); + Forward3D::Matrix3x3 invJacobMat; + calcInverseOfJacobianMatrix( JacobMat, detJacob, invJacobMat ); + const double dzdx = JacobMat.mat13 / JacobMat.mat11; + const double factorX = sqrt( 1.0 + dzdx * dzdx ); + const double dzdy = JacobMat.mat23 / JacobMat.mat22; + const double factorY = sqrt( 1.0 + dzdy * dzdy ); + + const double xCoordDifferenceOfSegment = (endCoord.X - startCoord.X) * lengthX * 0.5; + const double yCoordDifferenceOfSegment = (endCoord.Y - startCoord.Y) * lengthY * 0.5; + + if( integralXCompFirst ){ + //voltageDifference -= std::complex( xCoordDifferenceOfSegment * factorX, 0.0 ) * calcValueElectricFieldTangentialXFromAllEdges( elemID, 4, startCoord.X, startCoord.Y, -1.0 );// Since electric field is constant on edges + //voltageDifference -= std::complex( yCoordDifferenceOfSegment * factorY, 0.0 ) * calcValueElectricFieldTangentialYFromAllEdges( elemID, 4, endCoord.X, endCoord.Y, -1.0 );// Since electric field is constant on edges + calcInterpolatorVectorOfElectricFieldTangentialXFromAllEdges( elemID, 4, startCoord.X, startCoord.Y, -1.0, irhs, std::complex( - xCoordDifferenceOfSegment * factorX, 0.0 ) ); + calcInterpolatorVectorOfElectricFieldTangentialYFromAllEdges( elemID, 4, endCoord.X, endCoord.Y, -1.0, irhs, std::complex( - yCoordDifferenceOfSegment * factorY, 0.0 ) ); + }else{ + //voltageDifference -= std::complex( yCoordDifferenceOfSegment * factorY, 0.0 ) * calcValueElectricFieldTangentialYFromAllEdges( elemID, 4, startCoord.X, startCoord.Y, -1.0 );// Since electric field is constant on edges + //voltageDifference -= std::complex( xCoordDifferenceOfSegment * factorX, 0.0 ) * calcValueElectricFieldTangentialXFromAllEdges( elemID, 4, endCoord.X, endCoord.Y, -1.0 );// Since electric field is constant on edges + calcInterpolatorVectorOfElectricFieldTangentialYFromAllEdges( elemID, 4, startCoord.X, startCoord.Y, -1.0, irhs, std::complex( - yCoordDifferenceOfSegment * factorY, 0.0 ) ); + calcInterpolatorVectorOfElectricFieldTangentialXFromAllEdges( elemID, 4, endCoord.X, endCoord.Y, -1.0, irhs, std::complex( - xCoordDifferenceOfSegment * factorX, 0.0 ) ); + } + + const double areaFraction = 0.25 * fabs( (endCoord.X - startCoord.X) * (endCoord.Y - startCoord.Y) ); + const double area = m_MeshDataNonConformingHexaElement.calcAreaOfFace(elemID, 4); + double areaWithSign = 0.5 * area * areaFraction; + if( !rotationDirectionPlus ){ + areaWithSign *= -1.0; + } + //voltageDifference += calcValueRotatedElectricFieldNormal(elemID, 0.0, 0.0) * areaWithSign;// Since magnetic field is constant in the area + calcInterpolatorVectorOfRotatedElectricFieldNormal( elemID, 0.0, 0.0, irhs, areaWithSign ); + } } // Calculate interpolator vector of difference of voltage @@ -1091,8 +1199,65 @@ const MeshDataNonConformingHexaElement* Forward3DNonConformingHexaElement0thOrde // Calculate difference of voltage for brick element std::complex Forward3DNonConformingHexaElement0thOrder::calcVoltageDifference( const int nElem, const int* elememtsIncludingDipole, const CommonParameters::locationXY* localCoordinateValuesStartPoint, const CommonParameters::locationXY* localCoordinateValuesEndPoint ) const{ - OutputFiles::m_logFile << "Error : " << __FUNCTION__ << " is not implemented" << std::endl; - exit(1); + + std::complex voltageDifference = std::complex(0.0, 0.0); + for( int ielem = 0; ielem < nElem; ++ielem ){ + const int elemID = elememtsIncludingDipole[ielem]; + const CommonParameters::locationXY startCoord = localCoordinateValuesStartPoint[ielem]; + const CommonParameters::locationXY endCoord = localCoordinateValuesEndPoint[ielem]; + + bool rotationDirectionPlus = true; + CommonParameters::locationXY sharedPoint = {0.0, 0.0}; + const bool integralXCompFirst = doesIntegralXCompFirst( startCoord, endCoord, rotationDirectionPlus, sharedPoint ); + + const double lengthX = m_MeshDataNonConformingHexaElement.getEdgeLengthX(elemID); + const double lengthY = m_MeshDataNonConformingHexaElement.getEdgeLengthY(elemID); + Forward3D::Matrix3x3 JacobMat; + const double detJacob = calcJacobianMatrix( elemID, 0.0, 0.0, -1.0, JacobMat ); + Forward3D::Matrix3x3 invJacobMat; + calcInverseOfJacobianMatrix( JacobMat, detJacob, invJacobMat ); + const double dzdx = JacobMat.mat13 / JacobMat.mat11; + const double factorX = sqrt( 1.0 + dzdx * dzdx ); + const double dzdy = JacobMat.mat23 / JacobMat.mat22; + const double factorY = sqrt( 1.0 + dzdy * dzdy ); + + const double xCoordDifferenceOfSegment = (endCoord.X - startCoord.X) * lengthX * 0.5; + const double yCoordDifferenceOfSegment = (endCoord.Y - startCoord.Y) * lengthY * 0.5; + + if( integralXCompFirst ){ + voltageDifference -= std::complex( xCoordDifferenceOfSegment * factorX, 0.0 ) * calcValueElectricFieldTangentialXFromAllEdges( elemID, 4, startCoord.X, startCoord.Y, -1.0 );// Since electric field is constant on edges + voltageDifference -= std::complex( yCoordDifferenceOfSegment * factorY, 0.0 ) * calcValueElectricFieldTangentialYFromAllEdges( elemID, 4, endCoord.X, endCoord.Y, -1.0 );// Since electric field is constant on edges + }else{ + voltageDifference -= std::complex( yCoordDifferenceOfSegment * factorY, 0.0 ) * calcValueElectricFieldTangentialYFromAllEdges( elemID, 4, startCoord.X, startCoord.Y, -1.0 );// Since electric field is constant on edges + voltageDifference -= std::complex( xCoordDifferenceOfSegment * factorX, 0.0 ) * calcValueElectricFieldTangentialXFromAllEdges( elemID, 4, endCoord.X, endCoord.Y, -1.0 );// Since electric field is constant on edges + } + +#ifdef _DEBUG_WRITE + std::cout << "voltageDifference : " << voltageDifference << std::endl; +#endif + + const double areaFraction = 0.25 * fabs( (endCoord.X - startCoord.X) * (endCoord.Y - startCoord.Y) ); + const double area = m_MeshDataNonConformingHexaElement.calcAreaOfFace(elemID, 4); + double areaWithSign = 0.5 * area * areaFraction; + if( !rotationDirectionPlus ){ + areaWithSign *= -1.0; + } + +#ifdef _DEBUG_WRITE + const std::complex tmp1 = calcValueRotatedElectricFieldNormal(elemID, 0.0, 0.0); + const std::complex tmp2 = calcValueRotatedElectricFieldZDirection(elemID, 0.0, 0.0, -1.0); +#endif + + voltageDifference += calcValueRotatedElectricFieldNormal(elemID, 0.0, 0.0) * areaWithSign;// Since magnetic field is constant in the area + +#ifdef _DEBUG_WRITE + std::cout << "areaWithSign calcValueRotatedElectricFieldZDirection : " << areaWithSign << " " << tmp1 << std::endl; + std::cout << "voltageDifference total : " << voltageDifference << std::endl; +#endif + } + + return voltageDifference; + } // Calculate difference of voltage for tetra element @@ -1355,6 +1520,36 @@ void Forward3DNonConformingHexaElement0thOrder::calcArrayConvertIDGlobal2NonZero } +// Calculate 2D jacobian matrix for the Earth's surface +double Forward3DNonConformingHexaElement0thOrder::calc2DJacobianMatrixForEarthSurface( const int elemID, const double xi, const double eta, + Forward3D::Matrix2x2& jacobMat ) const{ + + // Array of reference coord xi values for each node + const double xiAtNode[4] = { -1.0, 1.0, 1.0, -1.0 }; + // Array of reference coord eta values for each node + const double etaAtNode[4] = { -1.0, -1.0, 1.0, 1.0 }; + + jacobMat.mat11 = 0.0; + jacobMat.mat12 = 0.0; + jacobMat.mat21 = 0.0; + jacobMat.mat22 = 0.0; + for( int i = 0; i < 4; ++i ){ + const double xiNode = xiAtNode[i]; + const double etaNode = etaAtNode[i]; + const double tmp1 = 0.25 * xiNode * (1.0 + etaNode * eta); + const double tmp2 = 0.25 * etaNode * (1.0 + xiNode * xi); + const int node = m_MeshDataNonConformingHexaElement.getNodesOfElements( elemID, i ); + const double xCoord = m_MeshDataNonConformingHexaElement.getXCoordinatesOfNodes(node); + const double yCoord = m_MeshDataNonConformingHexaElement.getYCoordinatesOfNodes(node); + jacobMat.mat11 += tmp1 * xCoord; + jacobMat.mat12 += tmp1 * yCoord; + jacobMat.mat21 += tmp2 * xCoord; + jacobMat.mat22 += tmp2 * yCoord; + } + + return jacobMat.mat11 * jacobMat.mat22 - jacobMat.mat12 * jacobMat.mat21; + +} // Make map converting master dofs after degeneration and MPC factors from slave dof after degeneration void Forward3DNonConformingHexaElement0thOrder::makeMapSlaveDofToMasterDofAndFactors(){ @@ -1667,6 +1862,62 @@ void Forward3DNonConformingHexaElement0thOrder::calcMPCConstants(){ } +// Add master dof and factor pair to m_slaveDofToMasterDofAndFactors +bool Forward3DNonConformingHexaElement0thOrder::doesIntegralXCompFirst( const CommonParameters::locationXY& startPoint, const CommonParameters::locationXY& endPoint, + bool& rotationDirectionPlus, CommonParameters::locationXY& sharedPoint ) const{ + + const double eps = 1.0e-12; + + bool intersectTwoEdges = true; + if( ( fabs( startPoint.X - 1.0 ) < eps && fabs( endPoint.Y - 1.0 ) < eps ) || + ( fabs( startPoint.Y - 1.0 ) < eps && fabs( endPoint.X - 1.0 ) < eps ) ){ + sharedPoint.X = 1.0; + sharedPoint.Y = 1.0; + }else if( ( fabs( startPoint.X - 1.0 ) < eps && fabs( endPoint.Y + 1.0 ) < eps ) || + ( fabs( startPoint.Y + 1.0 ) < eps && fabs( endPoint.X - 1.0 ) < eps ) ){ + sharedPoint.X = 1.0; + sharedPoint.Y = -1.0; + }else if( ( fabs( startPoint.X + 1.0 ) < eps && fabs( endPoint.Y + 1.0 ) < eps ) || + ( fabs( startPoint.Y + 1.0 ) < eps && fabs( endPoint.X + 1.0 ) < eps ) ){ + sharedPoint.X = -1.0; + sharedPoint.Y = -1.0; + }else if( ( fabs( startPoint.X + 1.0 ) < eps && fabs( endPoint.Y - 1.0 ) < eps ) || + ( fabs( startPoint.Y - 1.0 ) < eps && fabs( endPoint.X + 1.0 ) < eps ) ){ + sharedPoint.X = -1.0; + sharedPoint.Y = 1.0; + }else{// Not intersect two edges => integral X component first + sharedPoint.X = startPoint.X; + sharedPoint.Y = endPoint.Y; + intersectTwoEdges = false; + } + + const double outerProduct = ( sharedPoint.X - startPoint.X )*( endPoint.Y - startPoint.Y ) - ( sharedPoint.Y - startPoint.Y )*( endPoint.X - startPoint.X ); + if( outerProduct > 0 ){ + rotationDirectionPlus = true; + }else{ + rotationDirectionPlus = false; + } + + if( intersectTwoEdges ){ + if( sharedPoint.X * sharedPoint.Y < 0.0 ){ + if( rotationDirectionPlus ){ + return true; + }else{ + return false; + } + }else{ + if( rotationDirectionPlus ){ + return false; + }else{ + return true; + } + } + }else{// Not intersect two edges => integral X component first + return false; + } + +} + // Add master dof and factor pair to m_slaveDofToMasterDofAndFactors void Forward3DNonConformingHexaElement0thOrder::addMasterDofAndFactorPair( const int slaveDof, const int masterDof, const double factor ){ @@ -1910,6 +2161,31 @@ double Forward3DNonConformingHexaElement0thOrder::getShapeFuncRotatedZ( const do return tmp1 + tmp2; } +// Calculate jacobian matrix of the elements +double Forward3DNonConformingHexaElement0thOrder::get2DShapeFuncRotatedForEarthSurface( const double xi, const double eta, const int num, const Forward3D::Matrix2x2& invJacobMat ) const{ + + // Array of reference coord xi values for each edge + const double xiAtEdge[4] = { -9999.999, -9999.999, -1.0, 1.0 }; + // Array of reference coord eta values for each edge + const double etaAtEdge[4] = { -1.0, 1.0, -9999.999, -9999.999 }; + + switch( num ){ + case 0:// go through + case 1: + return 0.25 * etaAtEdge[num] * (invJacobMat.mat21*invJacobMat.mat12 - invJacobMat.mat11*invJacobMat.mat22 ); + break; + case 2:// go through + case 3: + return 0.25 * xiAtEdge[num] * (invJacobMat.mat22*invJacobMat.mat11 - invJacobMat.mat12*invJacobMat.mat21 ); + break; + default: + OutputFiles::m_logFile << "Error : Wrong number in " << __FUNCTION__ <<" : num = " << num << std::endl; + exit(1); + break; + } + +} + // Calculate jacobian matrix of the elements double Forward3DNonConformingHexaElement0thOrder::calcJacobianMatrix( const int elemID, const double xi, const double eta, const double zeta, Forward3D::Matrix3x3& JacobMat ) const{ diff --git a/src/Forward3DNonConformingHexaElement0thOrder.h b/src/Forward3DNonConformingHexaElement0thOrder.h index 03fd0cd..f7e5e2e 100644 --- a/src/Forward3DNonConformingHexaElement0thOrder.h +++ b/src/Forward3DNonConformingHexaElement0thOrder.h @@ -54,6 +54,9 @@ public: // Calculate Z component of rotated electric field virtual std::complex calcValueRotatedElectricFieldZDirection( const int iElem, const double xLocal, const double yLocal, const double zLocal ) const; + // Calculate X component of electric field only from the edges on the Earth's surface + std::complex calcValueRotatedElectricFieldNormal( const int iElem, const double xLocal, const double yLocal ) const; + // Calculate X component of electric field only from the edges on the Earth's surface virtual std::complex calcValueElectricFieldXDirectionFromEdgesOnEarthSurface( const int iElem, const int iFace, const double uCoord, const double vCoord ) const; @@ -93,6 +96,9 @@ public: // Calculate interpolator vector of Z component of rotated electric field virtual void calcInterpolatorVectorOfRotatedElectricFieldZDirection( const int iElem, const double xLocal, const double yLocal, const double zLocal, const int irhs, const std::complex& factor = std::complex(1.0,0.0) ); + // Calculate interpolator vector of X component of electric field only from the edges on the Earth's surface + void calcInterpolatorVectorOfRotatedElectricFieldNormal( const int iElem, const double xLocal, const double yLocal, const int irhs, const std::complex& factor = std::complex(1.0,0.0) ); + // Calculate interpolator vector of X component of electric field only from the edges on the Earth's surface virtual void calcInterpolatorVectorOfElectricFieldXDirectionFromEdgesOnEarthSurface( const int iElem, const int iFace, const double uCoord, const double vCoord, const int irhs, const std::complex& factor = std::complex(1.0,0.0) ); @@ -241,12 +247,20 @@ private: // Calculate array converting global edge IDs non-zero electric field values specified to the edges void calcArrayConvertIDGlobal2NonZeroValues(); + // Make map converting master dofs after degeneration and MPC factors from slave dof after degeneration + double calc2DJacobianMatrixForEarthSurface( const int elemID, const double xi, const double eta, + Forward3D::Matrix2x2& JacobMat ) const; + // Make map converting master dofs after degeneration and MPC factors from slave dof after degeneration void makeMapSlaveDofToMasterDofAndFactors(); // Calculate MPC constants void calcMPCConstants(); + // Add master dof and factor pair to m_slaveDofToMasterDofAndFactors + bool doesIntegralXCompFirst( const CommonParameters::locationXY& startPoint, const CommonParameters::locationXY& endPoint, + bool& rotationDirectionPlus, CommonParameters::locationXY& sharedPoint ) const; + // Add master dof and factor pair to m_slaveDofToMasterDofAndFactors void addMasterDofAndFactorPair( const int slaveDof, const int masterDof, const double factor ); @@ -268,6 +282,9 @@ private: // Get z component of shape function rotated for 0th order edge-based elements double getShapeFuncRotatedZ( const double xi, const double eta, const double zeta, const int num, const Forward3D::Matrix3x3& invJacobMat ) const; + // Calculate jacobian matrix of the elements + double get2DShapeFuncRotatedForEarthSurface( const double xi, const double eta, const int num, const Forward3D::Matrix2x2& invJacobMat ) const; + // Calculate jacobian matrix of the elements double calcJacobianMatrix( const int elemID, const double xi, const double eta, const double zeta, Forward3D::Matrix3x3& JacobMat ) const; diff --git a/src/MeshData.cpp b/src/MeshData.cpp index 416705f..7c4dffc 100644 --- a/src/MeshData.cpp +++ b/src/MeshData.cpp @@ -33,6 +33,7 @@ #include "CommonParameters.h" #include "ResistivityBlock.h" #include "OutputFiles.h" +#include "Util.h" // Constructer MeshData::MeshData(): @@ -347,18 +348,179 @@ double MeshData::calcDistance( const CommonParameters::locationXY& point0, cons } //// Decide whether specified elements share same edges -//bool MeshData::shareSameEdges( const int elemID1, const int elemID2 ) const{ -// -// if( elemID1 < 0 || elemID1 >= m_numElemTotal ){ -// OutputFiles::m_logFile << " Error : elemID1 is out of range in shareSameNodes. elemID1 = " << elemID1 << std::endl; -// exit(1); -// } -// -// if( elemID2 < 0 || elemID2 >= m_numElemTotal ){ -// OutputFiles::m_logFile << " Error : elemID2 is out of range in shareSameNodes. elemID2 = " << elemID2 << std::endl; -// exit(1); -// } -// -// -// -//} +bool MeshData::does1stSegmentContain2ndSegment( const CommonParameters::locationXY& startPointOf1stSegment, const CommonParameters::locationXY& endPointOf1stSegment, + const CommonParameters::locationXY& startPointOf2ndSegment, const CommonParameters::locationXY& endPointOf2ndSegment ) const{ + + const double dx1st = endPointOf1stSegment.X - startPointOf1stSegment.X; + const double dy1st = endPointOf1stSegment.Y - startPointOf1stSegment.Y; + const double outerProduct1 = dx1st * ( startPointOf2ndSegment.Y - startPointOf1stSegment.Y ) - dy1st * ( startPointOf2ndSegment.X - startPointOf1stSegment.X ); + const double outerProduct2 = dx1st * ( endPointOf2ndSegment.Y - startPointOf1stSegment.Y ) - dy1st * ( endPointOf2ndSegment.X - startPointOf1stSegment.X ); + + const double eps = 1.0e-9; + if( fabs(outerProduct1) < eps && fabs(outerProduct2) < eps ){ + // Two segmenta are located along a same line + const double innerProduct1 = ( startPointOf2ndSegment.X - startPointOf1stSegment.X ) * ( endPointOf1stSegment.X - startPointOf1stSegment.X ) + + ( startPointOf2ndSegment.Y - startPointOf1stSegment.Y ) * ( endPointOf1stSegment.Y - startPointOf1stSegment.Y ); + const double innerProduct2 = ( endPointOf2ndSegment.X - startPointOf1stSegment.X ) * ( endPointOf1stSegment.X - startPointOf1stSegment.X ) + + ( endPointOf2ndSegment.Y - startPointOf1stSegment.Y ) * ( endPointOf1stSegment.Y - startPointOf1stSegment.Y ); + const double denominator = pow(endPointOf1stSegment.X - startPointOf1stSegment.X, 2) + pow(endPointOf1stSegment.Y - startPointOf1stSegment.Y, 2); + if( innerProduct1 / denominator > - eps && innerProduct1 / denominator < 1.0 + eps && + innerProduct2 / denominator > - eps && innerProduct2 / denominator < 1.0 + eps ){ + return true; + } + } + + return false; +} + +// Function determine if two segments intersect or not +bool MeshData::intersectTwoSegments( const CommonParameters::locationXY& startPointOf1stSegment, const CommonParameters::locationXY& endPointOf1stSegment, + const CommonParameters::locationXY& startPointOf2ndSegment, const CommonParameters::locationXY& endPointOf2ndSegment ) const{ + + const double val1 = ( endPointOf1stSegment.Y - startPointOf1stSegment.Y ) * ( startPointOf2ndSegment.X - startPointOf1stSegment.X ) + ( endPointOf1stSegment.X - startPointOf1stSegment.X ) * ( startPointOf1stSegment.Y - startPointOf2ndSegment.Y ); + const double val2 = ( endPointOf1stSegment.Y - startPointOf1stSegment.Y ) * ( endPointOf2ndSegment.X - startPointOf1stSegment.X ) + ( endPointOf1stSegment.X - startPointOf1stSegment.X ) * ( startPointOf1stSegment.Y - endPointOf2ndSegment.Y ); + + const double eps = 1.0e-9; + + if( calcDistance(startPointOf1stSegment, startPointOf2ndSegment) < eps || + calcDistance(endPointOf1stSegment, startPointOf2ndSegment) < eps || + calcDistance(startPointOf1stSegment, endPointOf2ndSegment) < eps || + calcDistance(endPointOf1stSegment, endPointOf2ndSegment) < eps ){ + // Two segment share a node + return true; + } + + if( val1*val2 < eps ){ + const double val3 = ( endPointOf2ndSegment.Y - startPointOf2ndSegment.Y ) * ( startPointOf1stSegment.X - startPointOf2ndSegment.X ) + ( endPointOf2ndSegment.X - startPointOf2ndSegment.X ) * ( startPointOf2ndSegment.Y - startPointOf1stSegment.Y ); + const double val4 = ( endPointOf2ndSegment.Y - startPointOf2ndSegment.Y ) * ( endPointOf1stSegment.X - startPointOf2ndSegment.X ) + ( endPointOf2ndSegment.X - startPointOf2ndSegment.X ) * ( startPointOf2ndSegment.Y - endPointOf1stSegment.Y ); + if( fabs(val1*val2) < eps && fabs(val3*val4) < eps ){ + return false; + }else if( val3*val4 < eps ){ + return true; + } + //if( val3*val4 < eps ){ + // return true; + //} + } + + return false; + +} + +// Function determine if two lines overlap or not +bool MeshData::overlapTwoLines( const CommonParameters::locationXY& coord1stLine1, const CommonParameters::locationXY& coord1stLine2, + const CommonParameters::locationXY& coord2ndLine1, const CommonParameters::locationXY& coord2ndLine2 ) const{ + + const double eps = 1.0e-12; + if( fabs( coord1stLine2.X - coord1stLine1.X ) < eps ){// If the first line is parallel to Y direction + if( fabs( coord2ndLine2.X - coord2ndLine1.X ) < eps ){// If the second line is also parallel to Y direction + if( fabs( coord1stLine1.X - coord2ndLine1.X ) < eps && fabs( coord1stLine2.X - coord2ndLine2.X ) < eps ){ + return true; + }else{ + return false; + } + }else{// If the second line is not parallel to Y direction + return false; + } + } + + if( fabs( coord1stLine2.Y - coord1stLine1.Y ) < eps ){// If the first line is parallel to X direction + if( fabs( coord2ndLine2.Y - coord2ndLine1.Y ) < eps ){// If the second line is also parallel to X direction + if( fabs( coord1stLine1.Y - coord2ndLine1.Y ) < eps && fabs( coord1stLine2.Y - coord2ndLine2.Y ) < eps ){ + return true; + }else{ + return false; + } + }else{// If the second line is not parallel to X direction + return false; + } + } + + const double val1 = ( coord1stLine2.Y - coord1stLine1.Y )*( coord2ndLine2.X - coord2ndLine1.X ) - ( coord1stLine2.X - coord1stLine1.X )*( coord2ndLine2.Y - coord2ndLine1.Y ); + + const double val2 = ( coord1stLine2.Y - coord1stLine1.Y )*( coord2ndLine2.X - coord2ndLine1.X )*coord1stLine1.X + - ( coord1stLine2.X - coord1stLine1.X )*( coord2ndLine2.Y - coord2ndLine1.Y )*coord2ndLine1.X + + ( coord1stLine2.X - coord1stLine1.X )*( coord2ndLine2.X - coord2ndLine1.X )*( coord2ndLine1.Y - coord1stLine1.Y ); + + if( fabs(val1) < eps && fabs(val2) < eps ){ + return true; + } + + return false; + +} + +// Function determine if two segments overlap or not +bool MeshData::overlapTwoSegments( const CommonParameters::locationXY& startPointOf1stSegment, const CommonParameters::locationXY& endPointOf1stSegment, + const CommonParameters::locationXY& startPointOf2ndSegment, const CommonParameters::locationXY& endPointOf2ndSegment ) const{ + + if( !overlapTwoLines( startPointOf1stSegment, endPointOf1stSegment, startPointOf2ndSegment, endPointOf2ndSegment ) ){ + // Two lines don't overlap + return false; + } + + const double innerProduct1 = calcInnerProduct2D( startPointOf1stSegment, endPointOf1stSegment, startPointOf1stSegment, endPointOf1stSegment ); + const double innerProduct2 = calcInnerProduct2D( startPointOf1stSegment, endPointOf1stSegment, startPointOf1stSegment, startPointOf2ndSegment ); + const double innerProduct3 = calcInnerProduct2D( startPointOf1stSegment, endPointOf1stSegment, startPointOf1stSegment, endPointOf2ndSegment ); + + if( ( innerProduct2 < 0.0 && innerProduct3 < 0.0 ) || ( innerProduct2 > innerProduct1 && innerProduct3 > innerProduct1 ) ){ + return false; + } + + return true; + +} + +// Calculate inner product of two vectors +double MeshData::calcInnerProduct2D( const CommonParameters::locationXY& startCoordOf1stVec, const CommonParameters::locationXY& endCoordOf1stVec, + const CommonParameters::locationXY& startCoordOf2ndVec, const CommonParameters::locationXY& endCoordOf2ndVec) const{ + + CommonParameters::Vector3D vec1 = { endCoordOf1stVec.X - startCoordOf1stVec.X, endCoordOf1stVec.Y - startCoordOf1stVec.Y, 0.0 }; + CommonParameters::Vector3D vec2 = { endCoordOf2ndVec.X - startCoordOf2ndVec.X, endCoordOf2ndVec.Y - startCoordOf2ndVec.Y, 0.0 }; + + return calcInnerProduct( vec1, vec2 ); + +} + +// Calculate coordinates of intersection point of two lines +void MeshData::calcCoordOfIntersectionPointOfTwoLines( const CommonParameters::locationXY& coord1stLine1, const CommonParameters::locationXY& coord1stLine2, + const CommonParameters::locationXY& coord2ndLine1, const CommonParameters::locationXY& coord2ndLine2, CommonParameters::locationXY& coordIntersectionPoint) const{ + + const double eps = 1.0e-9; + if( calcDistance(coord1stLine1, coord2ndLine1) < eps || calcDistance(coord1stLine1, coord2ndLine2) < eps ){ + coordIntersectionPoint = coord1stLine1; + return; + } + if( calcDistance(coord1stLine2, coord2ndLine1) < eps || calcDistance(coord1stLine2, coord2ndLine2) < eps ){ + coordIntersectionPoint = coord1stLine2; + return; + } + + const double temp1 = ( coord1stLine2.Y - coord1stLine1.Y )*( coord2ndLine2.X - coord2ndLine1.X ) - ( coord1stLine2.X - coord1stLine1.X )*( coord2ndLine2.Y - coord2ndLine1.Y ); + if( fabs( temp1 ) < eps ){ + OutputFiles::m_logFile << " Error : Divide by zero in calculating X coordinate of intersection point of two lines !!" << std::endl; + exit(1); + } + + coordIntersectionPoint.X = ( coord1stLine2.Y - coord1stLine1.Y )*( coord2ndLine2.X - coord2ndLine1.X )*coord1stLine1.X + - ( coord1stLine2.X - coord1stLine1.X )*( coord2ndLine2.Y - coord2ndLine1.Y )*coord2ndLine1.X + + ( coord1stLine2.X - coord1stLine1.X )*( coord2ndLine2.X - coord2ndLine1.X )*( coord2ndLine1.Y - coord1stLine1.Y ); + coordIntersectionPoint.X /= temp1; + + const double temp2 = coord1stLine2.X - coord1stLine1.X; + const double temp3 = coord2ndLine2.X - coord2ndLine1.X; + + if( fabs( temp2 ) < 1.0e-8 && fabs( temp3 ) < 1.0e-8 ){ + OutputFiles::m_logFile << " Error : Divide by zero in calculating Y coordinate of intersection point of two lines !!" << std::endl; + exit(1); + } + + if( fabs( temp2 ) > fabs( temp3 ) ){ + coordIntersectionPoint.Y = ( coord1stLine2.Y - coord1stLine1.Y )/temp2*( coordIntersectionPoint.X - coord1stLine1.X ) + coord1stLine1.Y; + }else{ + coordIntersectionPoint.Y = ( coord2ndLine2.Y - coord2ndLine1.Y )/temp3*( coordIntersectionPoint.X - coord2ndLine1.X ) + coord2ndLine1.Y; + } + + return; + +} diff --git a/src/MeshData.h b/src/MeshData.h index 8623aed..0dacaf4 100644 --- a/src/MeshData.h +++ b/src/MeshData.h @@ -201,6 +201,30 @@ protected: // Calculate distanceof two points double calcDistance( const CommonParameters::locationXY& point0, const CommonParameters::locationXY& point1 ) const; + // Function determine if the 1st segment contains the 2nd segment + bool does1stSegmentContain2ndSegment( const CommonParameters::locationXY& startPointOf1stSegment, const CommonParameters::locationXY& endPointOf1stSegment, + const CommonParameters::locationXY& startPointOf2ndSegment, const CommonParameters::locationXY& endPointOf2ndSegment ) const; + + // Function determine if two segments intersect or not + bool intersectTwoSegments( const CommonParameters::locationXY& startPointOf1stSegment, const CommonParameters::locationXY& endPointOf1stSegment, + const CommonParameters::locationXY& startPointOf2ndSegment, const CommonParameters::locationXY& endPointOf2ndSegment ) const; + + // Function determine if two lines overlap or not + bool overlapTwoLines( const CommonParameters::locationXY& coord1stLine1, const CommonParameters::locationXY& coord1stLine2, + const CommonParameters::locationXY& coord2ndLine1, const CommonParameters::locationXY& coord2ndLine2 ) const; + + // Function determine if two segments overlap or not + bool overlapTwoSegments( const CommonParameters::locationXY& startPointOf1stSegment, const CommonParameters::locationXY& endPointOf1stSegment, + const CommonParameters::locationXY& startPointOf2ndSegment, const CommonParameters::locationXY& endPointOf2ndSegment ) const; + + // Calculate inner product of two vectors + double calcInnerProduct2D( const CommonParameters::locationXY& startCoordOf1stVec, const CommonParameters::locationXY& endCoordOf1stVec, + const CommonParameters::locationXY& startCoordOf2ndVec, const CommonParameters::locationXY& endCoordOf2ndVec) const; + + // Calculate coordinates of intersection point of two lines + void calcCoordOfIntersectionPointOfTwoLines( const CommonParameters::locationXY& coord1stLine1, const CommonParameters::locationXY& coord1stLine2, + const CommonParameters::locationXY& coord2ndLine1, const CommonParameters::locationXY& coord2ndLine2, CommonParameters::locationXY& coordIntersectionPoint) const; + }; #endif diff --git a/src/MeshDataNonConformingHexaElement.cpp b/src/MeshDataNonConformingHexaElement.cpp index 340481d..21de055 100644 --- a/src/MeshDataNonConformingHexaElement.cpp +++ b/src/MeshDataNonConformingHexaElement.cpp @@ -514,6 +514,279 @@ int MeshDataNonConformingHexaElement::findElementIncludingPointOnSurface( const } +// Find element including a point on the Y-Z plane and return element ID of 2D mesh +void MeshDataNonConformingHexaElement::findElementsIncludingDipoleOnSurface( const double locXStart, const double locYStart, const double locXEnd, const double locYEnd, + std::vector& elements, std::vector& localCoordXStartPoint, std::vector& localCoordYStartPoint, + std::vector& localCoordXEndPoint, std::vector& localCoordYEndPoint ) const{ + + const double thresholdVal = 1.0E-6; + + const CommonParameters::locationXY nodeCoordDipoleStart = { locXStart, locYStart }; + const CommonParameters::locationXY nodeCoordDipoleEnd = { locXEnd, locYEnd }; + + const double dipoleLength = hypot( ( locXEnd - locXStart ), ( locYEnd - locYStart ) ); + if( dipoleLength < thresholdVal ){ + OutputFiles::m_logFile << " Error : Length of dipole ( " << dipoleLength << " [m] ) is too small !! " << std::endl; + exit(1); + } + + double dummy(0.0); + int faceID(0); + + double localCoordStart[3] = { 0.0, 0.0, 0.0 }; + const int iElemStart = findElementIncludingPointOnSurface( locXStart, locYStart, faceID, + localCoordStart[0], localCoordStart[1], localCoordStart[2], false, false, dummy, dummy ); + assert(faceID == 4); + + double localCoordEnd[3] = { 0.0, 0.0, 0.0 }; + const int iElemEnd = findElementIncludingPointOnSurface( locXEnd, locYEnd, faceID, + localCoordEnd[0], localCoordEnd[1], localCoordEnd[2], false, false, dummy, dummy ); + assert(faceID == 4); + + std::vector intersectPointsAreadyFound[2]; + for( int iElem = 0; iElem < m_numElemOnLandSurface; ++iElem ){ + const int elemID = m_elemOnLandSurface[iElem]; + const int faceID = m_faceLandSurface[iElem]; + + const int nodeID0 = getNodesOfElements( elemID, m_faceID2NodeID[faceID][0] ); + const int nodeID1 = getNodesOfElements( elemID, m_faceID2NodeID[faceID][1] ); + const int nodeID2 = getNodesOfElements( elemID, m_faceID2NodeID[faceID][2] ); + const int nodeID3 = getNodesOfElements( elemID, m_faceID2NodeID[faceID][3] ); + + const CommonParameters::locationXY nodeCoord0 = { getXCoordinatesOfNodes( nodeID0 ), getYCoordinatesOfNodes( nodeID0 ) }; + const CommonParameters::locationXY nodeCoord1 = { getXCoordinatesOfNodes( nodeID1 ), getYCoordinatesOfNodes( nodeID1 ) }; + const CommonParameters::locationXY nodeCoord2 = { getXCoordinatesOfNodes( nodeID2 ), getYCoordinatesOfNodes( nodeID2 ) }; + const CommonParameters::locationXY nodeCoord3 = { getXCoordinatesOfNodes( nodeID3 ), getYCoordinatesOfNodes( nodeID3 ) }; + + CommonParameters::locationXY intersectPoints[4]; + if( overlapTwoSegments( nodeCoord0, nodeCoord1, nodeCoordDipoleStart, nodeCoordDipoleEnd ) ){// node0 - node1 + const double innerProduct1 = calcInnerProduct2D( nodeCoordDipoleStart, nodeCoordDipoleEnd, nodeCoordDipoleStart, nodeCoord0 ) / dipoleLength; + const double innerProduct2 = calcInnerProduct2D( nodeCoordDipoleStart, nodeCoordDipoleEnd, nodeCoordDipoleStart, nodeCoord1 ) / dipoleLength; + if( innerProduct1 < innerProduct2 ){ + if( innerProduct1 < 0.0 ){ + intersectPoints[0] = nodeCoordDipoleStart; + }else{ + intersectPoints[0] = nodeCoord0; + } + if( innerProduct2 > dipoleLength ){ + intersectPoints[1] = nodeCoordDipoleEnd; + }else{ + intersectPoints[1] = nodeCoord1; + } + }else{ + if( innerProduct2 < 0.0 ){ + intersectPoints[0] = nodeCoordDipoleStart; + }else{ + intersectPoints[0] = nodeCoord1; + } + if( innerProduct1 > dipoleLength ){ + intersectPoints[1] = nodeCoordDipoleEnd; + }else{ + intersectPoints[1] = nodeCoord0; + } + } + }else if( overlapTwoSegments( nodeCoord1, nodeCoord2, nodeCoordDipoleStart, nodeCoordDipoleEnd ) ){// node1 - node2 + const double innerProduct1 = calcInnerProduct2D( nodeCoordDipoleStart, nodeCoordDipoleEnd, nodeCoordDipoleStart, nodeCoord1 ) / dipoleLength; + const double innerProduct2 = calcInnerProduct2D( nodeCoordDipoleStart, nodeCoordDipoleEnd, nodeCoordDipoleStart, nodeCoord2 ) / dipoleLength; + if( innerProduct1 < innerProduct2 ){ + if( innerProduct1 < 0.0 ){ + intersectPoints[0] = nodeCoordDipoleStart; + }else{ + intersectPoints[0] = nodeCoord1; + } + if( innerProduct2 > dipoleLength ){ + intersectPoints[1] = nodeCoordDipoleEnd; + }else{ + intersectPoints[1] = nodeCoord2; + } + }else{ + if( innerProduct2 < 0.0 ){ + intersectPoints[0] = nodeCoordDipoleStart; + }else{ + intersectPoints[0] = nodeCoord2; + } + if( innerProduct1 > dipoleLength ){ + intersectPoints[1] = nodeCoordDipoleEnd; + }else{ + intersectPoints[1] = nodeCoord1; + } + } + }else if( overlapTwoSegments( nodeCoord2, nodeCoord3, nodeCoordDipoleStart, nodeCoordDipoleEnd ) ){// node2 - node3 + const double innerProduct1 = calcInnerProduct2D( nodeCoordDipoleStart, nodeCoordDipoleEnd, nodeCoordDipoleStart, nodeCoord2 ) / dipoleLength; + const double innerProduct2 = calcInnerProduct2D( nodeCoordDipoleStart, nodeCoordDipoleEnd, nodeCoordDipoleStart, nodeCoord3 ) / dipoleLength; + if( innerProduct1 < innerProduct2 ){ + if( innerProduct1 < 0.0 ){ + intersectPoints[0] = nodeCoordDipoleStart; + }else{ + intersectPoints[0] = nodeCoord2; + } + if( innerProduct2 > dipoleLength ){ + intersectPoints[1] = nodeCoordDipoleEnd; + }else{ + intersectPoints[1] = nodeCoord3; + } + }else{ + if( innerProduct2 < 0.0 ){ + intersectPoints[0] = nodeCoordDipoleStart; + }else{ + intersectPoints[0] = nodeCoord3; + } + if( innerProduct1 > dipoleLength ){ + intersectPoints[1] = nodeCoordDipoleEnd; + }else{ + intersectPoints[1] = nodeCoord2; + } + } + }else if( overlapTwoSegments( nodeCoord3, nodeCoord0, nodeCoordDipoleStart, nodeCoordDipoleEnd ) ){// node3 - node0 + const double innerProduct1 = calcInnerProduct2D( nodeCoordDipoleStart, nodeCoordDipoleEnd, nodeCoordDipoleStart, nodeCoord3 ) / dipoleLength; + const double innerProduct2 = calcInnerProduct2D( nodeCoordDipoleStart, nodeCoordDipoleEnd, nodeCoordDipoleStart, nodeCoord0 ) / dipoleLength; + if( innerProduct1 < innerProduct2 ){ + if( innerProduct1 < 0.0 ){ + intersectPoints[0] = nodeCoordDipoleStart; + }else{ + intersectPoints[0] = nodeCoord3; + } + if( innerProduct2 > dipoleLength ){ + intersectPoints[1] = nodeCoordDipoleEnd; + }else{ + intersectPoints[1] = nodeCoord0; + } + }else{ + if( innerProduct2 < 0.0 ){ + intersectPoints[0] = nodeCoordDipoleStart; + }else{ + intersectPoints[0] = nodeCoord0; + } + if( innerProduct1 > dipoleLength ){ + intersectPoints[1] = nodeCoordDipoleEnd; + }else{ + intersectPoints[1] = nodeCoord3; + } + } + }else{ + int icount(0); + if( intersectTwoSegments( nodeCoord0, nodeCoord1, nodeCoordDipoleStart, nodeCoordDipoleEnd ) ){ + calcCoordOfIntersectionPointOfTwoLines( nodeCoord0, nodeCoord1, nodeCoordDipoleStart, nodeCoordDipoleEnd, intersectPoints[icount] ); + ++icount; + } + if( intersectTwoSegments( nodeCoord1, nodeCoord2, nodeCoordDipoleStart, nodeCoordDipoleEnd ) ){ + calcCoordOfIntersectionPointOfTwoLines( nodeCoord1, nodeCoord2, nodeCoordDipoleStart, nodeCoordDipoleEnd, intersectPoints[icount] ); + bool duplicated(false); + for( int i = 0; i < icount; ++i ){ + if( calcDistance(intersectPoints[icount], intersectPoints[i]) < 1.0e-9 ){ + duplicated = true; + } + } + if(!duplicated){ + ++icount; + } + } + if( intersectTwoSegments( nodeCoord2, nodeCoord3, nodeCoordDipoleStart, nodeCoordDipoleEnd ) ){ + calcCoordOfIntersectionPointOfTwoLines( nodeCoord2, nodeCoord3, nodeCoordDipoleStart, nodeCoordDipoleEnd, intersectPoints[icount] ); + bool duplicated(false); + for( int i = 0; i < icount; ++i ){ + if( calcDistance(intersectPoints[icount], intersectPoints[i]) < 1.0e-9 ){ + duplicated = true; + } + } + if(!duplicated){ + ++icount; + } + } + if( intersectTwoSegments( nodeCoord3, nodeCoord0, nodeCoordDipoleStart, nodeCoordDipoleEnd ) ){ + calcCoordOfIntersectionPointOfTwoLines( nodeCoord3, nodeCoord0, nodeCoordDipoleStart, nodeCoordDipoleEnd, intersectPoints[icount] ); + bool duplicated(false); + for( int i = 0; i < icount; ++i ){ + if( calcDistance(intersectPoints[icount], intersectPoints[i]) < 1.0e-9 ){ + duplicated = true; + } + } + if(!duplicated){ + ++icount; + } + } + if( elemID == iElemStart && elemID == iElemEnd ){ + intersectPoints[0] = nodeCoordDipoleStart; + intersectPoints[1] = nodeCoordDipoleEnd; + }else if( icount == 0 ){ + continue; + }else if( icount == 1 ){ + if( elemID == iElemStart ){ + intersectPoints[1] = nodeCoordDipoleStart; + }else if( elemID == iElemEnd ){ + intersectPoints[1] = nodeCoordDipoleEnd; + }else{ + continue; + } + }else if( icount > 2 ){ + OutputFiles::m_logFile << " Error : Number of intersection points is larger than two: " << icount << std::endl; + exit(1); + } + + } + const double innerProduct = calcInnerProduct2D( nodeCoordDipoleStart, nodeCoordDipoleEnd, intersectPoints[0], intersectPoints[1] ); + if( fabs(innerProduct) < thresholdVal ){ + // Segment is too short + continue; + } + if( innerProduct < 0 ){ + const CommonParameters::locationXY temp = intersectPoints[0]; + intersectPoints[0] = intersectPoints[1]; + intersectPoints[1] = temp; + } + bool alreadyFound = false; + const int numSegments = static_cast( intersectPointsAreadyFound[0].size() ); + for( int iSeg = 0; iSeg < numSegments; ++iSeg ){ + const CommonParameters::locationXY coordStart = intersectPointsAreadyFound[0][iSeg]; + const CommonParameters::locationXY coordEnd = intersectPointsAreadyFound[1][iSeg]; + if( does1stSegmentContain2ndSegment(coordStart, coordEnd, intersectPoints[0], intersectPoints[1]) ){ + alreadyFound = true; + } + if( does1stSegmentContain2ndSegment(intersectPoints[0], intersectPoints[1], coordStart, coordEnd) ){ + alreadyFound = true; + } + } + if( alreadyFound ){ + // Segment has already been found + continue; + } + CommonParameters::locationXY localCoord = {0.0, 0.0}; + calcHorizontalLocalCoordinates( elemID, intersectPoints[0].X, intersectPoints[0].Y, localCoord.X, localCoord.Y ); + localCoordXStartPoint.push_back(localCoord.X); + localCoordYStartPoint.push_back(localCoord.Y); + calcHorizontalLocalCoordinates( elemID, intersectPoints[1].X, intersectPoints[1].Y, localCoord.X, localCoord.Y ); + localCoordXEndPoint.push_back(localCoord.X); + localCoordYEndPoint.push_back(localCoord.Y); + elements.push_back( elemID ); + for( int i = 0; i < 2; ++i ){ + intersectPointsAreadyFound[i].push_back( intersectPoints[i] ); + } + } + + if( elements.empty() ){ + OutputFiles::m_logFile << " Error : Could not find element including dipole ( " << locXStart << ", " << locYStart << " ) => ( " << locXEnd << ", " << locYEnd << " )." << std::endl; + exit(1); + } + + // For check + double dipoleLengthAccumulated(0.0); + std::vector::iterator itr0End = intersectPointsAreadyFound[0].end(); + std::vector::iterator itr0 = intersectPointsAreadyFound[0].begin(); + std::vector::iterator itr1 = intersectPointsAreadyFound[1].begin(); + while( itr0 != itr0End ){ + dipoleLengthAccumulated += calcDistance( *itr0, *itr1 ); + ++itr0; + ++itr1; + } + + if( fabs(dipoleLength - dipoleLengthAccumulated) / fabs(dipoleLength) > 0.01 ){ + OutputFiles::m_logFile << " Warning : Accumulated dipole length (" << dipoleLengthAccumulated + << ") is significantly different from the dipole length of the horizontal plane (" << dipoleLength + << ")" << std::endl; + } + +} + // Find element including a point on the Y-Z plane and return element ID of 2D mesh int MeshDataNonConformingHexaElement::findElementIncludingPointOnYZPlaneAndReturnElemID2D( const int iPlane, const double locY, const double locZ, double& xi, double& eta ) const{ @@ -1050,6 +1323,26 @@ double MeshDataNonConformingHexaElement::calcEdgeLengthFromElementAndEdgeBoundar } +// Get face index of neighbor element +double MeshDataNonConformingHexaElement::getEdgeLengthX( const int iElem ) const{ + + const int node0 = getNodesOfElements( iElem, 0 ); + const int node2 = getNodesOfElements( iElem, 2 ); + + return caldDiffXOfTwoNodes( node0, node2 ); + +} + +// Get length of the edges parallel to Y coordinate +double MeshDataNonConformingHexaElement::getEdgeLengthY( const int iElem ) const{ + + const int node0 = getNodesOfElements( iElem, 0 ); + const int node2 = getNodesOfElements( iElem, 2 ); + + return caldDiffYOfTwoNodes( node0, node2 ); + +} + // Get face index of neighbor element int MeshDataNonConformingHexaElement::getFaceIndexOfNeighborElement( const int iFace ) const{ diff --git a/src/MeshDataNonConformingHexaElement.h b/src/MeshDataNonConformingHexaElement.h index 2b5b7c6..f1a99b1 100644 --- a/src/MeshDataNonConformingHexaElement.h +++ b/src/MeshDataNonConformingHexaElement.h @@ -48,6 +48,10 @@ public: int findElementIncludingPointOnSurface( const double locX, const double locY, int& faceID, double& xi, double& eta, double& zeta, const bool useUpperElem, const bool modLoc, double& locXMod, double& locYMod ) const; + // Find element including a point on the Y-Z plane and return element ID of 2D mesh + void findElementsIncludingDipoleOnSurface( const double locXStart, const double locYStart, const double locXEnd, const double locYEnd, + std::vector& elements, std::vector& localCoordXStartPoint, std::vector& localCoordYStartPoint, std::vector& localCoordXEndPoint, std::vector& localCoordYEndPoint ) const; + // Find element including a point on the Y-Z plane and return element ID of 2D mesh int findElementIncludingPointOnYZPlaneAndReturnElemID2D( const int iPlane, const double locY, const double locZ, double& xi, double& eta ) const; @@ -129,6 +133,12 @@ public: // Calculate length of edges of elements on boundary planes double calcEdgeLengthFromElementAndEdgeBoundaryPlanes( const int iPlane, const int iElem, const int iEdge ) const; + // Get face index of neighbor element + double getEdgeLengthX( const int iElem ) const; + + // Get length of the edges parallel to Y coordinate + double getEdgeLengthY( const int iElem ) const; + // Get face index of neighbor element int getFaceIndexOfNeighborElement( const int iFace ) const; @@ -202,7 +212,7 @@ private: // Check whether the specified point is located in the specified element bool isLocatedInTheElement( const double x, const double y, const double z, const int iElem ) const; - + // Calculate local coordinates void calcLocalCoordinates( const int iElem, const double x, const double y, const double z, double& xi, double& eta, double& zeta ) const; diff --git a/src/MeshDataTetraElement.cpp b/src/MeshDataTetraElement.cpp index 27973f1..65e8783 100644 --- a/src/MeshDataTetraElement.cpp +++ b/src/MeshDataTetraElement.cpp @@ -726,15 +726,12 @@ void MeshDataTetraElement::findElementsIncludingDipoleOnSurface( const double lo std::cout << "icount : " << icount << std::endl; #endif - if( icount == 0 ){ - - if( elemID == iElemStart && elemID == iElemEnd ){ - intersectPoints[0] = nodeCoordDipoleStart; - intersectPoints[1] = nodeCoordDipoleEnd; - }else{ - continue; - } - + if( elemID == iElemStart && elemID == iElemEnd ){ + intersectPoints[0] = nodeCoordDipoleStart; + intersectPoints[1] = nodeCoordDipoleEnd; + } + else if( icount == 0 ){ + continue; }else if( icount == 1 ){ if( elemID == iElemStart ){ @@ -864,6 +861,13 @@ void MeshDataTetraElement::findElementsIncludingDipoleOnSurface( const double lo ++itr0; ++itr1; } + + if( fabs(dipoleLength - dipoleLengthAccumulated) / fabs(dipoleLength) > 0.01 ){ + OutputFiles::m_logFile << " Warning : Accumulated dipole length (" << dipoleLengthAccumulated + << ") is significantly different from the dipole length of the horizontal plane (" << dipoleLength + << ")" << std::endl; + } + } // Decide whether specified elements share same edges @@ -1382,191 +1386,6 @@ MeshDataTetraElement& MeshDataTetraElement::operator=(const MeshDataTetraElement } // Function determine if two segments intersect or not -bool MeshDataTetraElement::intersectTwoSegments( const CommonParameters::locationXY& startPointOf1stSegment, const CommonParameters::locationXY& endPointOf1stSegment, - const CommonParameters::locationXY& startPointOf2ndSegment, const CommonParameters::locationXY& endPointOf2ndSegment ) const{ - - const double val1 = ( endPointOf1stSegment.Y - startPointOf1stSegment.Y ) * ( startPointOf2ndSegment.X - startPointOf1stSegment.X ) + ( endPointOf1stSegment.X - startPointOf1stSegment.X ) * ( startPointOf1stSegment.Y - startPointOf2ndSegment.Y ); - const double val2 = ( endPointOf1stSegment.Y - startPointOf1stSegment.Y ) * ( endPointOf2ndSegment.X - startPointOf1stSegment.X ) + ( endPointOf1stSegment.X - startPointOf1stSegment.X ) * ( startPointOf1stSegment.Y - endPointOf2ndSegment.Y ); - - //const double eps = 1.0e-9; - -//#ifdef _DEBUG_WRITE -// std::cout << "startPointOf1stSegment : " << startPointOf1stSegment.X << " " << startPointOf1stSegment.Y << std::endl; // For debug -// std::cout << "endPointOf1stSegment : " << endPointOf1stSegment.X << " " << endPointOf1stSegment.Y << std::endl; // For debug -// std::cout << "startPointOf2ndSegment : " << startPointOf2ndSegment.X << " " << startPointOf2ndSegment.Y << std::endl; // For debug -// std::cout << "endPointOf2ndSegment : " << endPointOf2ndSegment.X << " " << endPointOf2ndSegment.Y << std::endl; // For debug -// std::cout << "val1 : " << ( endPointOf1stSegment.Y - startPointOf1stSegment.Y ) * ( startPointOf2ndSegment.X - startPointOf1stSegment.X ) << " " << ( endPointOf1stSegment.X - startPointOf1stSegment.X ) * ( startPointOf1stSegment.Y - startPointOf2ndSegment.Y ) << " " << val1 << std::endl; // For debug -// std::cout << "val2 : " << ( endPointOf1stSegment.Y - startPointOf1stSegment.Y ) * ( endPointOf2ndSegment.X - startPointOf1stSegment.X ) << " " << ( endPointOf1stSegment.X - startPointOf1stSegment.X ) * ( startPointOf1stSegment.Y - endPointOf2ndSegment.Y ) << " " << val2 << std::endl; // For debug -//#endif - - if( val1*val2 <= 0.0 ){ - - const double val3 = ( endPointOf2ndSegment.Y - startPointOf2ndSegment.Y ) * ( startPointOf1stSegment.X - startPointOf2ndSegment.X ) + ( endPointOf2ndSegment.X - startPointOf2ndSegment.X ) * ( startPointOf2ndSegment.Y - startPointOf1stSegment.Y ); - const double val4 = ( endPointOf2ndSegment.Y - startPointOf2ndSegment.Y ) * ( endPointOf1stSegment.X - startPointOf2ndSegment.X ) + ( endPointOf2ndSegment.X - startPointOf2ndSegment.X ) * ( startPointOf2ndSegment.Y - endPointOf1stSegment.Y ); - -//#ifdef _DEBUG_WRITE -// std::cout << "val3 : " << ( endPointOf2ndSegment.Y - startPointOf2ndSegment.Y ) * ( startPointOf1stSegment.X - startPointOf2ndSegment.X ) << " " << ( endPointOf2ndSegment.X - startPointOf2ndSegment.X ) * ( startPointOf2ndSegment.Y - startPointOf1stSegment.Y ) << " " << val3 << std::endl; // For debug -// std::cout << "val4 : " << ( endPointOf2ndSegment.Y - startPointOf2ndSegment.Y ) * ( endPointOf1stSegment.X - startPointOf2ndSegment.X ) << " " << ( endPointOf2ndSegment.X - startPointOf2ndSegment.X ) * ( startPointOf2ndSegment.Y - endPointOf1stSegment.Y ) << " " << val4 << std::endl; // For debug -//#endif - - if( fabs(val1*val2) < m_eps && fabs(val3*val4) < m_eps ){ - return false; - }else if( val3*val4 <= 0.0 ){ - return true; - } - - //if( val3*val4 <= 0.0 ){ - //return true; - //} - - } - - return false; - -} - -// Function determine if two lines overlap or not -bool MeshDataTetraElement::overlapTwoLines( const CommonParameters::locationXY& coord1stLine1, const CommonParameters::locationXY& coord1stLine2, - const CommonParameters::locationXY& coord2ndLine1, const CommonParameters::locationXY& coord2ndLine2 ) const{ - -//#ifdef _DEBUG_WRITE -// std::cout << "coord1stLine1 : " << coord1stLine1.X << " " << coord1stLine1.Y << std::endl; // For debug -// std::cout << "coord1stLine2 : " << coord1stLine2.X << " " << coord1stLine2.Y << std::endl; // For debug -// std::cout << "coord2ndLine1 : " << coord2ndLine1.X << " " << coord2ndLine1.Y << std::endl; // For debug -// std::cout << "coord2ndLine2 : " << coord2ndLine2.X << " " << coord2ndLine2.Y << std::endl; // For debug -//#endif - - //const double eps = 1.0e-6; - - if( fabs( coord1stLine2.X - coord1stLine1.X ) < m_eps ){// If the first line is parallel to Y direction - if( fabs( coord2ndLine2.X - coord2ndLine1.X ) < m_eps ){// If the second line is also parallel to Y direction - if( fabs( coord1stLine1.X - coord2ndLine1.X ) < m_eps && fabs( coord1stLine2.X - coord2ndLine2.X ) < m_eps ){ - return true; - }else{ - return false; - } - }else{// If the second line is not parallel to Y direction - return false; - } - } - - if( fabs( coord1stLine2.Y - coord1stLine1.Y ) < m_eps ){// If the first line is parallel to X direction - if( fabs( coord2ndLine2.Y - coord2ndLine1.Y ) < m_eps ){// If the second line is also parallel to X direction - if( fabs( coord1stLine1.Y - coord2ndLine1.Y ) < m_eps && fabs( coord1stLine2.Y - coord2ndLine2.Y ) < m_eps ){ - return true; - }else{ - return false; - } - }else{// If the second line is not parallel to X direction - return false; - } - } - - const double val1 = ( coord1stLine2.Y - coord1stLine1.Y )*( coord2ndLine2.X - coord2ndLine1.X ) - ( coord1stLine2.X - coord1stLine1.X )*( coord2ndLine2.Y - coord2ndLine1.Y ); - - const double val2 = ( coord1stLine2.Y - coord1stLine1.Y )*( coord2ndLine2.X - coord2ndLine1.X )*coord1stLine1.X - - ( coord1stLine2.X - coord1stLine1.X )*( coord2ndLine2.Y - coord2ndLine1.Y )*coord2ndLine1.X - + ( coord1stLine2.X - coord1stLine1.X )*( coord2ndLine2.X - coord2ndLine1.X )*( coord2ndLine1.Y - coord1stLine1.Y ); - -//#ifdef _DEBUG_WRITE -// std::cout << "1 : " << ( coord1stLine2.Y - coord1stLine1.Y )*( coord2ndLine2.X - coord2ndLine1.X )*coord1stLine1.X << std::endl; // For debug -// std::cout << "2 : " << - ( coord1stLine2.X - coord1stLine1.X )*( coord2ndLine2.Y - coord2ndLine1.Y )*coord2ndLine1.X << std::endl; // For debug -// std::cout << "3 : " << ( coord1stLine2.X - coord1stLine1.X )*( coord2ndLine2.X - coord2ndLine1.X )*( coord2ndLine1.Y - coord1stLine1.Y ) << std::endl; // For debug -// std::cout << "val1 : " << val1 << std::endl; // For debug -// std::cout << "val2 : " << val2 << std::endl; // For debug -//#endif - - if( fabs(val1) < m_eps && fabs(val2) < m_eps ){ - return true; - } - - return false; - -} - -// Function determine if two segments overlap or not -bool MeshDataTetraElement::overlapTwoSegments( const CommonParameters::locationXY& startPointOf1stSegment, const CommonParameters::locationXY& endPointOf1stSegment, - const CommonParameters::locationXY& startPointOf2ndSegment, const CommonParameters::locationXY& endPointOf2ndSegment ) const{ - - if( !overlapTwoLines( startPointOf1stSegment, endPointOf1stSegment, startPointOf2ndSegment, endPointOf2ndSegment ) ){ - // Two lines don't overlap - return false; - } - - const double innerProduct1 = calcInnerProduct2D( startPointOf1stSegment, endPointOf1stSegment, startPointOf1stSegment, endPointOf1stSegment ); - const double innerProduct2 = calcInnerProduct2D( startPointOf1stSegment, endPointOf1stSegment, startPointOf1stSegment, startPointOf2ndSegment ); - const double innerProduct3 = calcInnerProduct2D( startPointOf1stSegment, endPointOf1stSegment, startPointOf1stSegment, endPointOf2ndSegment ); - -//#ifdef _DEBUG_WRITE -// std::cout << "innerProduct1 : " << innerProduct1 << std::endl; // For debug -// std::cout << "innerProduct2 : " << innerProduct2 << std::endl; // For debug -// std::cout << "innerProduct3 : " << innerProduct3 << std::endl; // For debug -//#endif -// - if( ( innerProduct2 < 0.0 && innerProduct3 < 0.0 ) || ( innerProduct2 > innerProduct1 && innerProduct3 > innerProduct1 ) ){ - return false; - } - - return true; - -} - -// Calculate inner product of two vectors -double MeshDataTetraElement::calcInnerProduct2D( const CommonParameters::locationXY& startCoordOf1stVec, const CommonParameters::locationXY& endCoordOf1stVec, - const CommonParameters::locationXY& startCoordOf2ndVec, const CommonParameters::locationXY& endCoordOf2ndVec) const{ - - //return ( endCoordOf1stVec.X - startCoordOf1stVec.X )*( endCoordOf2ndVec.X - startCoordOf2ndVec.X )+( endCoordOf1stVec.Y - startCoordOf1stVec.Y )*( endCoordOf2ndVec.Y - startCoordOf2ndVec.Y ); - - CommonParameters::Vector3D vec1 = { endCoordOf1stVec.X - startCoordOf1stVec.X, endCoordOf1stVec.Y - startCoordOf1stVec.Y, 0.0 }; - CommonParameters::Vector3D vec2 = { endCoordOf2ndVec.X - startCoordOf2ndVec.X, endCoordOf2ndVec.Y - startCoordOf2ndVec.Y, 0.0 }; - - return calcInnerProduct( vec1, vec2 ); - -} - -// Calculate coordinates of intersection point of two lines -void MeshDataTetraElement::calcCoordOfIntersectionPointOfTwoLines( const CommonParameters::locationXY& coord1stLine1, const CommonParameters::locationXY& coord1stLine2, - const CommonParameters::locationXY& coord2ndLine1, const CommonParameters::locationXY& coord2ndLine2, CommonParameters::locationXY& coordIntersectionPoint) const{ - - const double temp1 = ( coord1stLine2.Y - coord1stLine1.Y )*( coord2ndLine2.X - coord2ndLine1.X ) - ( coord1stLine2.X - coord1stLine1.X )*( coord2ndLine2.Y - coord2ndLine1.Y ); - -//#ifdef _DEBUG_WRITE -// std::cout << "coord1stLine1 : " << coord1stLine1.X << " " << coord1stLine1.Y << std::endl; // For debug -// std::cout << "coord1stLine2 : " << coord1stLine2.X << " " << coord1stLine2.Y << std::endl; // For debug -// std::cout << "coord2ndLine1 : " << coord2ndLine1.X << " " << coord2ndLine1.Y << std::endl; // For debug -// std::cout << "coord2ndLine2 : " << coord2ndLine2.X << " " << coord2ndLine2.Y << std::endl; // For debug -// std::cout << "temp1 : " << ( coord1stLine2.Y - coord1stLine1.Y )*( coord2ndLine2.X - coord2ndLine1.X ) << " " << ( coord1stLine2.X - coord1stLine1.X )*( coord2ndLine2.Y - coord2ndLine1.Y ) << " " << temp1 << std::endl; // For debug -//#endif - - if( fabs( temp1 ) < 1.0e-12 ){ - OutputFiles::m_logFile << " Error : Divide by zero in calculating X coordinate of intersection point of two lines !!" << std::endl; - exit(1); - } - - coordIntersectionPoint.X = ( coord1stLine2.Y - coord1stLine1.Y )*( coord2ndLine2.X - coord2ndLine1.X )*coord1stLine1.X - - ( coord1stLine2.X - coord1stLine1.X )*( coord2ndLine2.Y - coord2ndLine1.Y )*coord2ndLine1.X - + ( coord1stLine2.X - coord1stLine1.X )*( coord2ndLine2.X - coord2ndLine1.X )*( coord2ndLine1.Y - coord1stLine1.Y ); - coordIntersectionPoint.X /= temp1; - - const double temp2 = coord1stLine2.X - coord1stLine1.X; - const double temp3 = coord2ndLine2.X - coord2ndLine1.X; - - if( fabs( temp2 ) < 1.0e-8 && fabs( temp3 ) < 1.0e-8 ){ - OutputFiles::m_logFile << " Error : Divide by zero in calculating Y coordinate of intersection point of two lines !!" << std::endl; - exit(1); - } - - if( fabs( temp2 ) > fabs( temp3 ) ){ - coordIntersectionPoint.Y = ( coord1stLine2.Y - coord1stLine1.Y )/temp2*( coordIntersectionPoint.X - coord1stLine1.X ) + coord1stLine1.Y; - }else{ - coordIntersectionPoint.Y = ( coord2ndLine2.Y - coord2ndLine1.Y )/temp3*( coordIntersectionPoint.X - coord2ndLine1.X ) + coord2ndLine1.Y; - } - - return; - -} - -// Function determine if then inputed point locate at the left of the surface of the lower element bool MeshDataTetraElement::locateLeftOfSegmentOnLandSurface( const CommonParameters::locationXY& point, const CommonParameters::locationXY& startPointOfSegment, const CommonParameters::locationXY& endPointOfSegment ) const{ diff --git a/src/MeshDataTetraElement.h b/src/MeshDataTetraElement.h index 3350aab..abd3f51 100644 --- a/src/MeshDataTetraElement.h +++ b/src/MeshDataTetraElement.h @@ -192,26 +192,6 @@ private: int m_edgeID2NodeID[6][2]; // Function determine if two segments intersect or not - bool intersectTwoSegments( const CommonParameters::locationXY& startPointOf1stSegment, const CommonParameters::locationXY& endPointOf1stSegment, - const CommonParameters::locationXY& startPointOf2ndSegment, const CommonParameters::locationXY& endPointOf2ndSegment ) const; - - // Function determine if two lines overlap or not - bool overlapTwoLines( const CommonParameters::locationXY& coord1stLine1, const CommonParameters::locationXY& coord1stLine2, - const CommonParameters::locationXY& coord2ndLine1, const CommonParameters::locationXY& coord2ndLine2 ) const; - - // Function determine if two segments overlap or not - bool overlapTwoSegments( const CommonParameters::locationXY& startPointOf1stSegment, const CommonParameters::locationXY& endPointOf1stSegment, - const CommonParameters::locationXY& startPointOf2ndSegment, const CommonParameters::locationXY& endPointOf2ndSegment ) const; - - // Calculate inner product of two vectors - double calcInnerProduct2D( const CommonParameters::locationXY& startCoordOf1stVec, const CommonParameters::locationXY& endCoordOf1stVec, - const CommonParameters::locationXY& startCoordOf2ndVec, const CommonParameters::locationXY& endCoordOf2ndVec) const; - - // Calculate coordinates of intersection point of two lines - void calcCoordOfIntersectionPointOfTwoLines( const CommonParameters::locationXY& coord1stLine1, const CommonParameters::locationXY& coord1stLine2, - const CommonParameters::locationXY& coord2ndLine1, const CommonParameters::locationXY& coord2ndLine2, CommonParameters::locationXY& coordIntersectionPoint) const; - - // Function determine if then inputed point locate at the left of the segment on the surface of the lower element bool locateLeftOfSegmentOnLandSurface( const CommonParameters::locationXY& point, const CommonParameters::locationXY& startPointOfSegment, const CommonParameters::locationXY& endPointOfSegment ) const; diff --git a/src/ObservedData.cpp b/src/ObservedData.cpp index 1a51d62..259b627 100644 --- a/src/ObservedData.cpp +++ b/src/ObservedData.cpp @@ -1436,6 +1436,8 @@ void ObservedData::outputCalculatedValuesOfAllStations() const{ } + fflush(OutputFiles::m_csvFile); + } // Calulate right-hand sides matrix consisting of interpolator vectors diff --git a/src/ObservedDataStationNMT.cpp b/src/ObservedDataStationNMT.cpp index 098a4aa..7aeeb35 100644 --- a/src/ObservedDataStationNMT.cpp +++ b/src/ObservedDataStationNMT.cpp @@ -253,7 +253,41 @@ void ObservedDataStationNMT::findElementsIncludingDipole(){ std::cout << m_elementsIncludingDipole[i] << " " << m_facesIncludingDipole[i] << " " << m_areaCoordinateValuesStartPoint[i].coord0 << " " << m_areaCoordinateValuesStartPoint[i].coord1 << " " << m_areaCoordinateValuesStartPoint[i].coord2 << std::endl; } #endif + }else if( ( AnalysisControl::getInstance() )->getTypeOfMesh() == MeshData::NONCONFORMING_HEXA ){// Non-conforming deformed hexahedral mesh + const MeshDataNonConformingHexaElement* const ptrMeshDataNonConformingHexaElement = ( AnalysisControl::getInstance() )->getPointerOfMeshDataNonConformingHexaElement(); + + std::vector localCoordXStartPoint; + std::vector localCoordYStartPoint; + std::vector localCoordXEndPoint; + std::vector localCoordYEndPoint; + std::vector elementsIncludingDipole; + + ptrMeshDataNonConformingHexaElement->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( 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 }else{// Hexa mesh const MeshDataBrickElement* const ptrMeshDataBrickElement = ( AnalysisControl::getInstance() )->getPointerOfMeshDataBrickElement(); @@ -593,6 +627,25 @@ double ObservedDataStationNMT::getZCoordOfPoint( const int num ) const{ } return ptrMeshDataTetraElement->calcZCoordOfPointOnFace( iElem, iFace, areaCoord ); + } + else if( ( AnalysisControl::getInstance() )->getTypeOfMesh() == MeshData::NONCONFORMING_HEXA ){ + + const MeshDataNonConformingHexaElement* const ptrMeshDataNonConformingHexaElement = ( AnalysisControl::getInstance() )->getPointerOfMeshDataNonConformingHexaElement(); + + int faceID(0); + double dummy(0.0); + double localCoordX(0.0); + double localCoordY(0.0); + double localCoordZ(0.0); + if( num == 0 ){ + iElem = ptrMeshDataNonConformingHexaElement->findElementIncludingPointOnSurface( + m_location.startPoint.X, m_location.startPoint.Y, faceID, localCoordX, localCoordY, localCoordZ, false, false, dummy, dummy ); + }else{ + iElem = ptrMeshDataNonConformingHexaElement->findElementIncludingPointOnSurface( + m_location.endPoint.X, m_location.endPoint.Y, faceID, localCoordX, localCoordY, localCoordZ, false, false, dummy, dummy ); + } + return ptrMeshDataNonConformingHexaElement->calcZCoordOfPointOnFace(iElem, 4, localCoordX, localCoordY); + } else{ diff --git a/src/ObservedDataStationNMT2.cpp b/src/ObservedDataStationNMT2.cpp index 454e578..b9cac80 100644 --- a/src/ObservedDataStationNMT2.cpp +++ b/src/ObservedDataStationNMT2.cpp @@ -357,7 +357,42 @@ void ObservedDataStationNMT2::findElementsIncludingDipoles(){ std::cout << m_elementsIncludingDipole[1][i] << " " << m_facesIncludingDipole[1][i] << " " << m_areaCoordinateValuesEndPoint[1][i].coord0 << " " << m_areaCoordinateValuesEndPoint[1][i].coord1 << " " << m_areaCoordinateValuesStartPoint[1][i].coord2 << std::endl; } #endif - + }else if( ( AnalysisControl::getInstance() )->getTypeOfMesh() == MeshData::NONCONFORMING_HEXA ){// Non-conforming deformed hexahedral mesh + const MeshDataNonConformingHexaElement* const ptrMeshDataNonConformingHexaElement = ( AnalysisControl::getInstance() )->getPointerOfMeshDataNonConformingHexaElement(); + for( int idipole = 0; idipole < 2; ++idipole ){ + std::vector localCoordXStartPoint; + std::vector localCoordYStartPoint; + std::vector localCoordXEndPoint; + std::vector localCoordYEndPoint; + std::vector elementsIncludingDipole; + ptrMeshDataNonConformingHexaElement->findElementsIncludingDipoleOnSurface( + m_location[idipole].startPoint.X, m_location[idipole].startPoint.Y, m_location[idipole].endPoint.X, m_location[idipole].endPoint.Y, + elementsIncludingDipole, localCoordXStartPoint, localCoordYStartPoint, localCoordXEndPoint, localCoordYEndPoint ); + m_numElementsIncludingDipole[idipole] = static_cast( elementsIncludingDipole.size() ); + if( m_numElementsIncludingDipole[idipole] > 0 ){ + m_elementsIncludingDipole[idipole] = new int[ m_numElementsIncludingDipole[idipole] ]; + m_localCoordinateValuesStartPoint[idipole] = new CommonParameters::locationXY[ m_numElementsIncludingDipole[idipole] ]; + m_localCoordinateValuesEndPoint[idipole] = new CommonParameters::locationXY[ m_numElementsIncludingDipole[idipole] ]; + } + for( int ielem = 0; ielem < m_numElementsIncludingDipole[idipole]; ++ielem ){ + m_elementsIncludingDipole[idipole][ielem] = elementsIncludingDipole[ielem]; + m_localCoordinateValuesStartPoint[idipole][ielem].X = localCoordXStartPoint[ielem]; + m_localCoordinateValuesStartPoint[idipole][ielem].Y = localCoordYStartPoint[ielem]; + m_localCoordinateValuesEndPoint[idipole][ielem].X = localCoordXEndPoint[ielem]; + m_localCoordinateValuesEndPoint[idipole][ielem].Y = localCoordYEndPoint[ielem]; + } + } +#ifdef _DEBUG_WRITE + std::cout << "m_stationID " << m_stationID << std::endl; + std::cout << "m_elementsIncludingDipole[0] m_localCoordinateValuesStartPoint[0].x m_localCoordinateValuesStartPoint[0].y m_localCoordinateValuesEndPoint[0].x m_localCoordinateValuesEndPoint[0].y" << std::endl; + for( int ielem = 0; ielem < m_numElementsIncludingDipole[0]; ++ielem ){ + std::cout << m_elementsIncludingDipole[0][ielem] << " " << m_localCoordinateValuesStartPoint[0][ielem].X << " " << m_localCoordinateValuesStartPoint[0][ielem].Y << " " << m_localCoordinateValuesEndPoint[0][ielem].X << " " << m_localCoordinateValuesEndPoint[0][ielem].Y << std::endl; + } + std::cout << "m_elementsIncludingDipole[1] m_localCoordinateValuesStartPoint[1].x m_localCoordinateValuesStartPoint[1].y m_localCoordinateValuesEndPoint[1].x m_localCoordinateValuesEndPoint[1].y" << std::endl; + for( int ielem = 0; ielem < m_numElementsIncludingDipole[1]; ++ielem ){ + std::cout << m_elementsIncludingDipole[1][ielem] << " " << m_localCoordinateValuesStartPoint[1][ielem].X << " " << m_localCoordinateValuesStartPoint[1][ielem].Y << " " << m_localCoordinateValuesEndPoint[1][ielem].X << " " << m_localCoordinateValuesEndPoint[1][ielem].Y << std::endl; + } +#endif }else{// Hexa mesh const MeshDataBrickElement* const ptrMeshDataBrickElement = ( AnalysisControl::getInstance() )->getPointerOfMeshDataBrickElement(); @@ -910,6 +945,25 @@ double ObservedDataStationNMT2::getZCoordOfPoint( const int iDipole , const int } return ptrMeshDataTetraElement->calcZCoordOfPointOnFace( iElem, iFace, areaCoord ); + } + else if( ( AnalysisControl::getInstance() )->getTypeOfMesh() == MeshData::NONCONFORMING_HEXA ){ + + const MeshDataNonConformingHexaElement* const ptrMeshDataNonConformingHexaElement = ( AnalysisControl::getInstance() )->getPointerOfMeshDataNonConformingHexaElement(); + + int faceID(0); + double dummy(0.0); + double localCoordX(0.0); + double localCoordY(0.0); + double localCoordZ(0.0); + if( num == 0 ){ + iElem = ptrMeshDataNonConformingHexaElement->findElementIncludingPointOnSurface( + m_location[iDipole].startPoint.X, m_location[iDipole].startPoint.Y, faceID, localCoordX, localCoordY, localCoordZ, false, false, dummy, dummy ); + }else{ + iElem = ptrMeshDataNonConformingHexaElement->findElementIncludingPointOnSurface( + m_location[iDipole].endPoint.X, m_location[iDipole].endPoint.Y, faceID, localCoordX, localCoordY, localCoordZ, false, false, dummy, dummy ); + } + return ptrMeshDataNonConformingHexaElement->calcZCoordOfPointOnFace(iElem, 4, localCoordX, localCoordY); + } else{