From 7097a4eaa0a32e0d02207521941157bda8968b05 Mon Sep 17 00:00:00 2001 From: Ruben Undheim Date: Mon, 13 Aug 2018 09:26:34 +0200 Subject: New upstream version 0.0.35+ds.1 --- AppCSXCAD/CMakeLists.txt | 10 +- CSXCAD/CMakeLists.txt | 6 +- CSXCAD/src/CMakeLists.txt | 1 + CSXCAD/src/CSPrimBox.cpp | 6 +- CSXCAD/src/CSPrimitives.cpp | 12 +- CSXCAD/src/CSPrimitives.h | 3 +- CSXCAD/src/CSTransform.cpp | 5 + CSXCAD/src/CSTransform.h | 3 + QCSXCAD/CMakeLists.txt | 4 +- QCSXCAD/VTKPrimitives.cpp | 18 +- openEMS/CMakeLists.txt | 17 +- openEMS/Common/processmodematch.cpp | 12 +- openEMS/FDTD/excitation.cpp | 55 +++-- openEMS/FDTD/extensions/engine_ext_excitation.cpp | 12 +- openEMS/FDTD/extensions/engine_ext_tfsf.cpp | 4 +- openEMS/FDTD/openems_fdtd_mpi.cpp | 7 +- openEMS/main.cpp | 8 +- openEMS/matlab/DelayFidelity.m | 93 +++++++++ openEMS/matlab/RunOpenEMS.m | 2 +- openEMS/matlab/RunOpenEMS_MPI.m | 2 +- openEMS/matlab/Tutorials/CylindricalWave_CC.m | 2 +- openEMS/matlab/Tutorials/RadarUWBTutorial.m | 236 ++++++++++++++++++++++ openEMS/matlab/Tutorials/Simple_Patch_Antenna.m | 53 ++--- openEMS/matlab/Tutorials/StripLine2MSL.m | 133 ++++++++++++ openEMS/matlab/h5readatt_octave.cc | 21 +- openEMS/matlab/plotRefl.m | 145 +++++++++++++ openEMS/matlab/setup.m | 20 +- openEMS/nf2ff/CMakeLists.txt | 3 + openEMS/openEMS_MPI.pro | 5 - openEMS/openems.cpp | 9 +- openEMS/tools/hdf5_file_reader.cpp | 4 +- openEMS/tools/hdf5_file_writer.cpp | 14 +- 32 files changed, 777 insertions(+), 148 deletions(-) create mode 100644 openEMS/matlab/DelayFidelity.m create mode 100644 openEMS/matlab/Tutorials/RadarUWBTutorial.m create mode 100644 openEMS/matlab/Tutorials/StripLine2MSL.m create mode 100644 openEMS/matlab/plotRefl.m delete mode 100644 openEMS/openEMS_MPI.pro diff --git a/AppCSXCAD/CMakeLists.txt b/AppCSXCAD/CMakeLists.txt index d5b6e5a..d5a72c9 100644 --- a/AppCSXCAD/CMakeLists.txt +++ b/AppCSXCAD/CMakeLists.txt @@ -6,7 +6,7 @@ ELSE() SET( CMAKE_BUILD_TYPE Release CACHE STRING "Set to either \"Release\" or \"Debug\"" ) ENDIF() -PROJECT( AppCSXCAD CXX) +PROJECT( AppCSXCAD CXX C) cmake_minimum_required(VERSION 2.8) @@ -19,7 +19,7 @@ IF(EXISTS ${PROJECT_SOURCE_DIR}/localConfig.cmake) ENDIF() # default -set(VERSION "v0.2.1") +set(VERSION "v0.2.2") # add git revision IF(EXISTS ${PROJECT_SOURCE_DIR}/.git ) @@ -95,6 +95,10 @@ find_path(QCSXCAD_INCLUDE_DIR message(STATUS "QCSXCAD_INCLUDE_DIR: ${QCSXCAD_INCLUDE_DIR}" ) INCLUDE_DIRECTORIES( ${QCSXCAD_INCLUDE_DIR} ) +# hdf5 +find_package(HDF5 1.8 COMPONENTS C HL REQUIRED) +INCLUDE_DIRECTORIES (${HDF5_INCLUDE_DIRS}) + # vtk if (WIN32) find_package(VTK 6.1 REQUIRED) @@ -149,6 +153,8 @@ endif() TARGET_LINK_LIBRARIES( AppCSXCAD ${CSXCAD_LIBRARIES} ${QCSXCAD_LIBRARIES} + ${HDF5_LIBRARIES} + ${HDF5_HL_LIBRARIES} ${QT_LIBRARIES} ${vtk_LIBS} ) diff --git a/CSXCAD/CMakeLists.txt b/CSXCAD/CMakeLists.txt index a349287..9d516d5 100644 --- a/CSXCAD/CMakeLists.txt +++ b/CSXCAD/CMakeLists.txt @@ -6,14 +6,14 @@ ELSE() SET( CMAKE_BUILD_TYPE Release CACHE STRING "Set to either \"Release\" or \"Debug\"" ) ENDIF() -PROJECT(CSXCAD CXX) +PROJECT(CSXCAD CXX C) cmake_minimum_required(VERSION 2.8) # default set(LIB_VERSION_MAJOR 0) set(LIB_VERSION_MINOR 6) -set(LIB_VERSION_PATCH 1) +set(LIB_VERSION_PATCH 2) set(LIB_VERSION_STRING ${LIB_VERSION_MAJOR}.${LIB_VERSION_MINOR}.${LIB_VERSION_PATCH}) set(VERSION "v${LIB_VERSION_STRING}") @@ -91,7 +91,7 @@ find_package(TinyXML REQUIRED) ADD_DEFINITIONS( -DTIXML_USE_STL ) find_package(HDF5 1.8 COMPONENTS C HL REQUIRED) -INCLUDE_DIRECTORIES (${HDF5_INCLUDE_DIR}) +INCLUDE_DIRECTORIES (${HDF5_INCLUDE_DIRS}) link_directories(${HDF5_LIBRARY_DIRS}) # hdf5 compat ADD_DEFINITIONS( -DH5_USE_16_API ) diff --git a/CSXCAD/src/CMakeLists.txt b/CSXCAD/src/CMakeLists.txt index 676546a..8903a60 100644 --- a/CSXCAD/src/CMakeLists.txt +++ b/CSXCAD/src/CMakeLists.txt @@ -89,6 +89,7 @@ TARGET_LINK_LIBRARIES( CSXCAD ${fparser_LIBRARIES} ${TinyXML_LIBRARIES} ${HDF5_LIBRARIES} + ${HDF5_HL_LIBRARIES} CGAL ${Boost_LIBRARIES} ${vtk_LIBS} diff --git a/CSXCAD/src/CSPrimBox.cpp b/CSXCAD/src/CSPrimBox.cpp index de31b31..e4acf77 100644 --- a/CSXCAD/src/CSPrimBox.cpp +++ b/CSXCAD/src/CSPrimBox.cpp @@ -56,9 +56,6 @@ CSPrimBox::~CSPrimBox() bool CSPrimBox::GetBoundBox(double dBoundBox[6], bool PreserveOrientation) { -// if ( (m_MeshType!=m_PrimCoordSystem) && (m_PrimCoordSystem!=UNDEFINED_CS)) -// std::cerr << "GetBoundBox::GetBoundBox: Warning: The bounding box for this object is not calculated properly... " << std::endl; - const double* start = m_Coords[0].GetCoords(m_MeshType); const double* stop = m_Coords[1].GetCoords(m_MeshType); @@ -81,6 +78,9 @@ bool CSPrimBox::GetBoundBox(double dBoundBox[6], bool PreserveOrientation) dBoundBox[2*i]=dBoundBox[2*i+1]; dBoundBox[2*i+1]=help; } + if ( (m_MeshType!=m_PrimCoordSystem) && (m_PrimCoordSystem!=UNDEFINED_CS)) + // if the box is defined in a coordinate system other than the expected one, this BB is invalid + return false; return true; } diff --git a/CSXCAD/src/CSPrimitives.cpp b/CSXCAD/src/CSPrimitives.cpp index e15d2cb..f852328 100644 --- a/CSXCAD/src/CSPrimitives.cpp +++ b/CSXCAD/src/CSPrimitives.cpp @@ -123,6 +123,13 @@ void CSPrimitives::Init() m_BoundBoxValid = false; } +CSTransform* CSPrimitives::GetTransform() +{ + if (m_Transform==NULL) + m_Transform = new CSTransform(clParaSet); + return m_Transform; +} + void CSPrimitives::SetProperty(CSProperties *prop) { if ((clProperty!=NULL) && (clProperty!=prop)) @@ -146,8 +153,9 @@ int CSPrimitives::IsInsideBox(const double *boundbox) return 0; // unable to decide with an invalid bounding box if ((this->GetBoundBoxCoordSystem()!=UNDEFINED_CS) && (this->GetBoundBoxCoordSystem()!=this->GetCoordInputType())) return 0; // unable to decide if coordinate system do not match - if (this->GetTransform()!=NULL) - return 0; // unable to decide if a transformation is used + if (m_Transform!=NULL) + if (m_Transform->HasTransform()) + return 0; // unable to decide if a transformation is used for (int i=0;i<3;++i) { diff --git a/CSXCAD/src/CSPrimitives.h b/CSXCAD/src/CSPrimitives.h index 4f68a96..545c152 100644 --- a/CSXCAD/src/CSPrimitives.h +++ b/CSXCAD/src/CSPrimitives.h @@ -177,7 +177,8 @@ public: //! Read the coordinate system for this primitive (may be different to the input mesh type) \sa GetCoordInputType CoordinateSystem GetCoordinateSystem() const {return m_PrimCoordSystem;} - CSTransform* GetTransform() const {return m_Transform;} + //! Get the CSTransform if it exists already or create a new one + CSTransform* GetTransform(); //! Show status of this primitve virtual void ShowPrimitiveStatus(std::ostream& stream); diff --git a/CSXCAD/src/CSTransform.cpp b/CSXCAD/src/CSTransform.cpp index 2e32481..d71cb4e 100644 --- a/CSXCAD/src/CSTransform.cpp +++ b/CSXCAD/src/CSTransform.cpp @@ -73,6 +73,11 @@ void CSTransform::Reset() MakeUnitMatrix(m_Inv_TMatrix); } +bool CSTransform::HasTransform() +{ + return (m_TransformList.size()>0); +} + void CSTransform::Invert() { //make sure the inverse matrix is up to date... diff --git a/CSXCAD/src/CSTransform.h b/CSXCAD/src/CSTransform.h index 8293219..8fcad90 100644 --- a/CSXCAD/src/CSTransform.h +++ b/CSXCAD/src/CSTransform.h @@ -85,6 +85,9 @@ public: void Reset(); + //! Check if this CSTransform has any transformations + bool HasTransform(); + //! All subsequent operations will be occur before the previous operations (not the default). void SetPreMultiply() {m_PostMultiply=false;} //! All subsequent operations will be after the previous operations (default). diff --git a/QCSXCAD/CMakeLists.txt b/QCSXCAD/CMakeLists.txt index 32c80cd..73b96a6 100644 --- a/QCSXCAD/CMakeLists.txt +++ b/QCSXCAD/CMakeLists.txt @@ -22,7 +22,7 @@ ENDIF() # default set(LIB_VERSION_MAJOR 0) set(LIB_VERSION_MINOR 6) -set(LIB_VERSION_PATCH 1) +set(LIB_VERSION_PATCH 2) set(LIB_VERSION_STRING ${LIB_VERSION_MAJOR}.${LIB_VERSION_MINOR}.${LIB_VERSION_PATCH}) set(VERSION "v${LIB_VERSION_STRING}") @@ -107,7 +107,7 @@ endif() message(STATUS "Found package VTK. Using version " ${VTK_VERSION}) include(${VTK_USE_FILE}) -INCLUDE_DIRECTORIES (${VTK_INCLUDE_DIR}) +INCLUDE_DIRECTORIES (${VTK_INCLUDE_DIRS}) # Qt SET(RESOURCES resources.qrc) diff --git a/QCSXCAD/VTKPrimitives.cpp b/QCSXCAD/VTKPrimitives.cpp index 5610a9e..61c7c20 100644 --- a/QCSXCAD/VTKPrimitives.cpp +++ b/QCSXCAD/VTKPrimitives.cpp @@ -778,7 +778,7 @@ vtkActor* VTKPrimitives::AddPolyData(vtkPolyData* polydata, double *dRGB, double filter->SetTransform(vtrans); #if VTK_MAJOR_VERSION>=6 - m_PolyDataCollection->AddInputData(filter->GetOutput()); + m_PolyDataCollection->AddInputConnection(filter->GetOutputPort()); #else m_PolyDataCollection->AddInput(filter->GetOutput()); #endif @@ -812,7 +812,7 @@ vtkActor* VTKPrimitives::AddPolyData(vtkAlgorithmOutput* polydata_port, double * filter->SetTransform(vtrans); #if VTK_MAJOR_VERSION>=6 - m_PolyDataCollection->AddInputData(filter->GetOutput()); + m_PolyDataCollection->AddInputConnection(filter->GetOutputPort()); #else m_PolyDataCollection->AddInput(filter->GetOutput()); #endif @@ -855,7 +855,7 @@ void VTKPrimitives::WritePolyData2File(const char* filename, double scale) if (scale==1.0) { #if VTK_MAJOR_VERSION>=6 - writer->SetInputData(m_PolyDataCollection->GetOutput()); + writer->SetInputConnection(m_PolyDataCollection->GetOutputPort()); #else writer->SetInput(m_PolyDataCollection->GetOutput()); #endif @@ -867,7 +867,7 @@ void VTKPrimitives::WritePolyData2File(const char* filename, double scale) vtkTransformPolyDataFilter *transformFilter = vtkTransformPolyDataFilter::New(); #if VTK_MAJOR_VERSION>=6 - transformFilter->SetInputData(m_PolyDataCollection->GetOutput()); + transformFilter->SetInputConnection(m_PolyDataCollection->GetOutputPort()); #else transformFilter->SetInput(m_PolyDataCollection->GetOutput()); #endif @@ -875,7 +875,7 @@ void VTKPrimitives::WritePolyData2File(const char* filename, double scale) transformFilter->SetTransform(transform); #if VTK_MAJOR_VERSION>=6 - writer->SetInputData(transformFilter->GetOutput()); + writer->SetInputConnection(transformFilter->GetOutputPort()); #else writer->SetInput(transformFilter->GetOutput()); #endif @@ -895,7 +895,7 @@ void VTKPrimitives::WritePolyData2STL(const char* filename, double scale) vtkTriangleFilter* filter = vtkTriangleFilter::New(); #if VTK_MAJOR_VERSION>=6 - filter->SetInputData(m_PolyDataCollection->GetOutput()); + filter->SetInputConnection(m_PolyDataCollection->GetOutputPort()); #else filter->SetInput(m_PolyDataCollection->GetOutput()); #endif @@ -917,7 +917,7 @@ void VTKPrimitives::WritePolyData2STL(const char* filename, double scale) transformFilter->SetTransform(transform); #if VTK_MAJOR_VERSION>=6 - writer->SetInputData(transformFilter->GetOutput()); + writer->SetInputConnection(transformFilter->GetOutputPort()); #else writer->SetInput(transformFilter->GetOutput()); #endif @@ -937,7 +937,7 @@ void VTKPrimitives::WritePolyData2PLY(const char* filename, double scale) vtkTriangleFilter* filter = vtkTriangleFilter::New(); #if VTK_MAJOR_VERSION>=6 - filter->SetInputData(m_PolyDataCollection->GetOutput()); + filter->SetInputConnection(m_PolyDataCollection->GetOutputPort()); #else filter->SetInput(m_PolyDataCollection->GetOutput()); #endif @@ -960,7 +960,7 @@ void VTKPrimitives::WritePolyData2PLY(const char* filename, double scale) transformFilter->SetTransform(transform); #if VTK_MAJOR_VERSION>=6 - writer->SetInputData(transformFilter->GetOutput()); + writer->SetInputConnection(transformFilter->GetOutputPort()); #else writer->SetInput(transformFilter->GetOutput()); #endif diff --git a/openEMS/CMakeLists.txt b/openEMS/CMakeLists.txt index 5434138..898f280 100644 --- a/openEMS/CMakeLists.txt +++ b/openEMS/CMakeLists.txt @@ -6,13 +6,13 @@ ELSE() SET( CMAKE_BUILD_TYPE Release CACHE STRING "Set to either \"Release\" or \"Debug\"" ) ENDIF() -PROJECT(openEMS CXX) +PROJECT(openEMS CXX C) cmake_minimum_required(VERSION 2.8) # default set(LIB_VERSION_MAJOR 0) set(LIB_VERSION_MINOR 0) -set(LIB_VERSION_PATCH 34) +set(LIB_VERSION_PATCH 35) set(LIB_VERSION_STRING ${LIB_VERSION_MAJOR}.${LIB_VERSION_MINOR}.${LIB_VERSION_PATCH}) set(VERSION "v${LIB_VERSION_STRING}") @@ -21,7 +21,7 @@ IF(EXISTS ${PROJECT_SOURCE_DIR}/localConfig.cmake) include(${PROJECT_SOURCE_DIR}/localConfig.cmake) ENDIF() -set(VERSION "v0.0.34") +set(VERSION "v0.0.35") # add git revision IF(EXISTS ${PROJECT_SOURCE_DIR}/.git ) @@ -109,11 +109,8 @@ ADD_DEFINITIONS( -DTIXML_USE_STL ) # hdf5 find_package(HDF5 1.8 COMPONENTS C HL REQUIRED) -INCLUDE_DIRECTORIES (${HDF5_INCLUDE_DIR}) -link_directories(${HDF5_LIBRARY_DIRS}) - -# hdf5 compat -ADD_DEFINITIONS( -DH5_USE_16_API ) +INCLUDE_DIRECTORIES (${HDF5_INCLUDE_DIRS}) +link_directories(${HDF5_LIBRARIES}) # boost find_package(Boost 1.46 COMPONENTS @@ -158,6 +155,8 @@ set(SOURCES openems.cpp ) +set(PUB_HEADERS openems.h) + # libs ADD_SUBDIRECTORY( tools ) ADD_SUBDIRECTORY( FDTD ) @@ -181,6 +180,7 @@ TARGET_LINK_LIBRARIES( openEMS ${fparser_LIBRARIES} tinyxml ${HDF5_LIBRARIES} + ${HDF5_HL_LIBRARIES} ${Boost_LIBRARIES} ${vtk_LIBS} ${MPI_LIBRARIES} @@ -204,5 +204,6 @@ if (UNIX) PERMISSIONS OWNER_WRITE OWNER_READ GROUP_READ WORLD_READ OWNER_EXECUTE GROUP_EXECUTE WORLD_EXECUTE) endif() endif () +INSTALL(FILES ${PUB_HEADERS} DESTINATION include/openEMS) INSTALL( DIRECTORY matlab DESTINATION share/openEMS ) # TODO mpi, tarball, debug, release diff --git a/openEMS/Common/processmodematch.cpp b/openEMS/Common/processmodematch.cpp index 9c71c7a..620257d 100644 --- a/openEMS/Common/processmodematch.cpp +++ b/openEMS/Common/processmodematch.cpp @@ -135,22 +135,23 @@ void ProcessModeMatch::InitProcess() m_ModeDist[n] = Create2DArray(m_numLines); } + bool dualMesh = m_ModeFieldType==1; unsigned int pos[3] = {0,0,0}; double discLine[3] = {0,0,0}; double gridDelta = 1; // 1 -> mode-matching function is definied in drawing units... double var[7]; pos[m_ny] = start[m_ny]; - discLine[m_ny] = Op->GetDiscLine(m_ny,pos[m_ny],m_dualMesh); + discLine[m_ny] = Op->GetDiscLine(m_ny,pos[m_ny],dualMesh); double norm = 0; double area = 0; for (unsigned int posP = 0; posPGetDiscLine(nP,pos[nP],m_dualMesh); + discLine[nP] = Op->GetDiscLine(nP,pos[nP],dualMesh); for (unsigned int posPP = 0; posPPGetDiscLine(nPP,pos[nPP],m_dualMesh); + discLine[nPP] = Op->GetDiscLine(nPP,pos[nPP],dualMesh); var[0] = discLine[0] * gridDelta; // x var[1] = discLine[1] * gridDelta; // y @@ -169,7 +170,7 @@ void ProcessModeMatch::InitProcess() var[5] = sqrt(pow(discLine[0],2)+pow(discLine[2],2)) * gridDelta; // r var[6] = asin(1)-atan(var[2]/var[3]); //theta (t) } - area = Op->GetNodeArea(m_ny,pos,m_dualMesh); + area = Op->GetNodeArea(m_ny,pos,dualMesh); for (int n=0; n<2; ++n) { m_ModeDist[n][posP][posPP] = m_ModeParser[n]->Eval(var); //calc mode template @@ -227,6 +228,7 @@ double* ProcessModeMatch::CalcMultipleIntegrals() double field = 0; double purity = 0; double area = 0; + bool dualMesh = m_ModeFieldType==1; int nP = (m_ny+1)%3; int nPP = (m_ny+2)%3; @@ -242,7 +244,7 @@ double* ProcessModeMatch::CalcMultipleIntegrals() for (unsigned int posPP = 0; posPPGetNodeArea(m_ny,pos,m_dualMesh); + area = Op->GetNodeArea(m_ny,pos,dualMesh); if (m_ModeFieldType==0) m_Eng_Interface->GetEField(pos,out); if (m_ModeFieldType==1) diff --git a/openEMS/FDTD/excitation.cpp b/openEMS/FDTD/excitation.cpp index f095b53..4e1ba85 100644 --- a/openEMS/FDTD/excitation.cpp +++ b/openEMS/FDTD/excitation.cpp @@ -137,7 +137,7 @@ unsigned int Excitation::GetMaxExcitationTimestep() const { FDTD_FLOAT maxAmp=0; unsigned int maxStep=0; - for (unsigned int n=1; nmaxAmp) { @@ -152,7 +152,7 @@ void Excitation::CalcGaussianPulsExcitation(double f0, double fc, int nTS) { if (dT==0) return; - Length = (unsigned int)(2.0 * 9.0/(2.0*PI*fc) / dT); + Length = (unsigned int)ceil(2.0 * 9.0/(2.0*PI*fc) / dT); if (Length>(unsigned int)nTS) { cerr << "Operator::CalcGaussianPulsExcitation: Requested excitation pusle would be " << Length << " timesteps or " << Length * dT << " s long. Cutting to max number of timesteps!" << endl; @@ -160,13 +160,11 @@ void Excitation::CalcGaussianPulsExcitation(double f0, double fc, int nTS) } delete[] Signal_volt; delete[] Signal_curr; - Signal_volt = new FDTD_FLOAT[Length+1]; - Signal_curr = new FDTD_FLOAT[Length+1]; - Signal_volt[0]=0.0; - Signal_curr[0]=0.0; - for (unsigned int n=1; nVolt_delay[n]; exc_pos *= (exc_pos>0); exc_pos %= p; - exc_pos *= (exc_pos<=(int)length); + exc_pos *= (exc_pos<(int)length); ny = m_Op_Exc->Volt_dir[n]; pos[0]=m_Op_Exc->Volt_index[0][n]; pos[1]=m_Op_Exc->Volt_index[1][n]; @@ -71,7 +71,7 @@ void Engine_Ext_Excitation::Apply2Voltages() exc_pos = numTS - (int)m_Op_Exc->Volt_delay[n]; exc_pos *= (exc_pos>0); exc_pos %= p; - exc_pos *= (exc_pos<=(int)length); + exc_pos *= (exc_pos<(int)length); ny = m_Op_Exc->Volt_dir[n]; pos[0]=m_Op_Exc->Volt_index[0][n]; pos[1]=m_Op_Exc->Volt_index[1][n]; @@ -87,7 +87,7 @@ void Engine_Ext_Excitation::Apply2Voltages() exc_pos = numTS - (int)m_Op_Exc->Volt_delay[n]; exc_pos *= (exc_pos>0); exc_pos %= p; - exc_pos *= (exc_pos<=(int)length); + exc_pos *= (exc_pos<(int)length); ny = m_Op_Exc->Volt_dir[n]; pos[0]=m_Op_Exc->Volt_index[0][n]; pos[1]=m_Op_Exc->Volt_index[1][n]; @@ -124,7 +124,7 @@ void Engine_Ext_Excitation::Apply2Current() exc_pos = numTS - (int)m_Op_Exc->Curr_delay[n]; exc_pos *= (exc_pos>0); exc_pos %= p; - exc_pos *= (exc_pos<=(int)length); + exc_pos *= (exc_pos<(int)length); ny = m_Op_Exc->Curr_dir[n]; pos[0]=m_Op_Exc->Curr_index[0][n]; pos[1]=m_Op_Exc->Curr_index[1][n]; @@ -141,7 +141,7 @@ void Engine_Ext_Excitation::Apply2Current() exc_pos = numTS - (int)m_Op_Exc->Curr_delay[n]; exc_pos *= (exc_pos>0); exc_pos %= p; - exc_pos *= (exc_pos<=(int)length); + exc_pos *= (exc_pos<(int)length); ny = m_Op_Exc->Curr_dir[n]; pos[0]=m_Op_Exc->Curr_index[0][n]; pos[1]=m_Op_Exc->Curr_index[1][n]; @@ -157,7 +157,7 @@ void Engine_Ext_Excitation::Apply2Current() exc_pos = numTS - (int)m_Op_Exc->Curr_delay[n]; exc_pos *= (exc_pos>0); exc_pos %= p; - exc_pos *= (exc_pos<=(int)length); + exc_pos *= (exc_pos<(int)length); ny = m_Op_Exc->Curr_dir[n]; pos[0]=m_Op_Exc->Curr_index[0][n]; pos[1]=m_Op_Exc->Curr_index[1][n]; diff --git a/openEMS/FDTD/extensions/engine_ext_tfsf.cpp b/openEMS/FDTD/extensions/engine_ext_tfsf.cpp index cd3b107..9942656 100644 --- a/openEMS/FDTD/extensions/engine_ext_tfsf.cpp +++ b/openEMS/FDTD/extensions/engine_ext_tfsf.cpp @@ -44,7 +44,7 @@ void Engine_Ext_TFSF::DoPostVoltageUpdates() { if ( numTS < n ) m_DelayLookup[n]=0; - else if ((numTS-n > length) && (p==0)) + else if ((numTS-n >= length) && (p==0)) m_DelayLookup[n]=0; else m_DelayLookup[n] = numTS - n; @@ -134,7 +134,7 @@ void Engine_Ext_TFSF::DoPostCurrentUpdates() { if ( numTS < n ) m_DelayLookup[n]=0; - else if ((numTS-n > length) && (p==0)) + else if ((numTS-n >= length) && (p==0)) m_DelayLookup[n]=0; else m_DelayLookup[n] = numTS - n; diff --git a/openEMS/FDTD/openems_fdtd_mpi.cpp b/openEMS/FDTD/openems_fdtd_mpi.cpp index eaab59a..964bdd4 100644 --- a/openEMS/FDTD/openems_fdtd_mpi.cpp +++ b/openEMS/FDTD/openems_fdtd_mpi.cpp @@ -210,7 +210,7 @@ bool openEMS_FDTD_MPI::SetupMPI() if (numProcs!=m_NumProc) { if (m_MyID==0) - cerr << "openEMS_FDTD_MPI::SetupMPI: Error: Requested splits require " << numProcs << " processes, but only " << m_NumProc << " were found! Exit! " << endl; + cerr << "openEMS_FDTD_MPI::SetupMPI: Error: Requested splits require " << numProcs << " processes, but " << m_NumProc << " were found! Exit! " << endl; exit(10); } @@ -256,8 +256,8 @@ bool openEMS_FDTD_MPI::SetupMPI() grid->AddDiscLine(2, m_Original_Grid->GetLine(2,n) ); m_MPI_Op->SetSplitPos(0,m_SplitNumber[0].at(i)); - m_MPI_Op->SetSplitPos(1,m_SplitNumber[1].at(i)); - m_MPI_Op->SetSplitPos(2,m_SplitNumber[2].at(i)); + m_MPI_Op->SetSplitPos(1,m_SplitNumber[1].at(j)); + m_MPI_Op->SetSplitPos(2,m_SplitNumber[2].at(k)); if (i>0) m_MPI_Op->SetNeighborDown(0,procTable[i-1][j][k]); @@ -398,6 +398,7 @@ bool openEMS_FDTD_MPI::SetupProcessing() { //type is integral processing --> disable! Needs to be fixed! cerr << "openEMS_FDTD_MPI::SetupProcessing(): Warning: Processing: " << proc->GetName() << " occures multiple times and is being deactivated..." << endl; + cerr << "openEMS_FDTD_MPI::SetupProcessing(): Note: Processing: Make sure that there are no splits inside probes or sources." << endl; deactivate = true; rename = false; } diff --git a/openEMS/main.cpp b/openEMS/main.cpp index 1d03481..8a8000f 100644 --- a/openEMS/main.cpp +++ b/openEMS/main.cpp @@ -68,8 +68,12 @@ int main(int argc, char *argv[]) } int EC = FDTD.ParseFDTDSetup(argv[1]); + if(!EC) { + cerr << "openEMS - ParseFDTDSetup failed." << endl; + exit(1); + } EC = FDTD.SetupFDTD(); - if (EC) return EC; + if (EC) exit(EC); FDTD.RunFDTD(); #ifdef MPI_SUPPORT @@ -77,5 +81,5 @@ int main(int argc, char *argv[]) MPI::Finalize(); #endif - return 0; + exit(0); } diff --git a/openEMS/matlab/DelayFidelity.m b/openEMS/matlab/DelayFidelity.m new file mode 100644 index 0000000..3cc0ae2 --- /dev/null +++ b/openEMS/matlab/DelayFidelity.m @@ -0,0 +1,93 @@ +function [delay, fidelity, nf2ff_out] = DelayFidelity(nf2ff, port, path, weight_theta, weight_phi, theta, phi, f_0, f_c, varargin) +% [delay, fidelity] = DelayFidelity(nf2ff, port, path, theta, phi, f_lo, f_hi, varargin) +% +% +% This function calculates the time delay from the source port to the phase center of the antenna and the fidelity. +% The fidelity is the similarity between the excitation pulse and the radiated pulse (normalized scalar product). +% The resolution of the delay will be equal to or better than ((f_0 + f_c)*Oversampling)^-1 when using Gaussian excitation. +% Oversampling is an input parameter to InitFDTD. The rows of delay and fidelity correspond to theta and the columns to phi. +% +% input: +% nf2ff: return value of CreateNF2FFBox. +% port: return value of AddLumpedPort +% path: path of the simulation results. +% weight_theta: weight if the E_theta component +% weight_phi: eight of the E_phi component +% -> with both (possibly complex) parameters any polarization can be examined +% theta: theta values to be simulated +% phi: phi values to be simulated +% f_0: center frequency of SetGaussExcite +% f_c: cutoff frequency of SetGaussExcite +% +% variable input: +% 'Center': phase center of the antenna for CalcNF2FF +% 'Radius': radius for CalcNF2FF +% 'Mode': mode CalcNF2FF +% +% example: +% theta = [-180:10:180] * pi / 180; +% phi = [0, 90] * pi / 180; +% [delay, fidelity] = DelayFidelity2(nf2ff, port, Sim_Path, sin(tilt), cos(tilt), theta, phi, f_0, f_c, 'Mode', 1); +% figure +% polar(theta.', delay(:,1) * 3e11); % delay in mm +% figure +% polar(theta', (fidelity(:,1)-0.95)/0.05); % last 5 percent of fidelity +% +% Author: Georg Michel + +C0 = 299792458; +center = [0, 0, 0]; +radius = 1; +nf2ff_mode = 0; + +for n=1:2:numel(varargin) + if (strcmp(varargin{n},'Center')==1); + center = varargin{n+1}; + elseif (strcmp(varargin{n},'Radius')==1); + radius = varargin{n+1}; + elseif (strcmp(varargin{n},'Mode')==1); + nf2ff_mode = varargin{n+1}; + end +end + + +port_ut = load(fullfile(path, port.U_filename)); +port_it = load(fullfile(path, port.I_filename)); +dt = port_ut(2,1) - port_ut(1,1); +fftsize = 2^(nextpow2(size(port_ut)(1)) + 1); +df = 1 / (dt * fftsize); +uport = fft(port_ut(:, 2), fftsize)(1:fftsize/2+1); +iport = fft(port_it(:, 2), fftsize)(1:fftsize/2+1); +fport = df * (0:fftsize/2); +f_ind = find(fport > (f_0 - f_c ) & fport < (f_0 + f_c)); +disp(["frequencies: ", num2str(numel(f_ind))]); +exc_f = uport.' + iport.' * port.Feed_R; %excitation in freq domain +exc_f(!f_ind) = 0; +exc_f /= sqrt(exc_f * exc_f'); % normalization (transposing also conjugates) + +nf2ff = CalcNF2FF(nf2ff, path, fport(f_ind), theta, phi, ... + 'Center', center, 'Radius', radius, 'Mode', nf2ff_mode); +radfield = weight_theta * cell2mat(nf2ff.E_theta) + weight_phi * cell2mat(nf2ff.E_phi); % rows: theta(f1), columns: phi(f1), phi(f2), ...phi(fn) +radfield = reshape(radfield, [length(nf2ff.theta), length(nf2ff.phi), length(nf2ff.freq)]); +correction = reshape(exp(-2i*pi*nf2ff.r/C0*nf2ff.freq), 1,1,numel(nf2ff.freq)); %dimensions: theta, phi, frequencies +radfield = radfield./correction; % correct for radius delay +% normalize radfield +radnorm = sqrt(dot(radfield, radfield, 3)); +radfield ./= radnorm; + +%initialize radiated field in fully populated frequency domain +rad_f = zeros([numel(nf2ff.theta), numel(nf2ff.phi), numel(fport)]); +rad_f(:, :, f_ind) = radfield; % assign selected frequencies +exc_f = reshape(exc_f, [1,1,numel(exc_f)]); %make exc_f confomant with rad_f + +cr_f = rad_f .* conj(exc_f); % calculate cross correlation +% calculate the cross correlation in time domain (analytic signal) +cr = ifft(cr_f(:, :, 1:end-1), [], 3) * (numel(fport) -1); % twice the FFT normalization (sqrt^2) because product of two normalized functions +%search for the maxiumum of the envelope +[fidelity, delay_ind] = max(abs(cr), [], 3); +delay = (delay_ind - 1) * dt * 2; % double time step because of single-sided FFT +nf2ff_out = nf2ff; %possibly needed for plotting the far field and other things +disp(["DelayFidelity: delay resolution = ", num2str(dt*2e9), "ns"]); +return; + + diff --git a/openEMS/matlab/RunOpenEMS.m b/openEMS/matlab/RunOpenEMS.m index a4947e3..3894adc 100644 --- a/openEMS/matlab/RunOpenEMS.m +++ b/openEMS/matlab/RunOpenEMS.m @@ -21,7 +21,7 @@ function RunOpenEMS(Sim_Path, Sim_File, opts, Settings) % --engine=fastest fastest available engine (default) % --engine=basic basic FDTD engine % --engine=sse engine using sse vector extensions -% --engine=sse_compressed engine using compressed operator + sse vector extensions +% --engine=sse-compressed engine using compressed operator + sse vector extensions % --engine=MPI engine using compressed operator + sse vector extensions + MPI parallel processing % --engine=multithreaded engine using compressed operator + sse vector extensions + MPI + multithreading % --numThreads= Force use n threads for multithreaded engine diff --git a/openEMS/matlab/RunOpenEMS_MPI.m b/openEMS/matlab/RunOpenEMS_MPI.m index 495f7e5..779b038 100644 --- a/openEMS/matlab/RunOpenEMS_MPI.m +++ b/openEMS/matlab/RunOpenEMS_MPI.m @@ -87,7 +87,7 @@ end if isfield(Settings.MPI,'Hosts') disp(['Running remote openEMS_MPI in working dir: ' work_path]); - [status] = system(['mpiexec -host ' HostList ' -n ' int2str(NrProc) ' -wdir ' work_path ' ' Settings.MPI.Binary ' ' Sim_File ' ' opts ' ' append_unix]); + [status] = system(['mpiexec ' Settings.MPI.GlobalArgs ' -host ' HostList ' -n ' int2str(NrProc) ' -wdir ' work_path ' ' Settings.MPI.Binary ' ' Sim_File ' ' opts ' ' append_unix]); else disp('Running local openEMS_MPI'); [status] = system(['mpiexec ' Settings.MPI.GlobalArgs ' -n ' int2str(NrProc) ' ' Settings.MPI.Binary ' ' Sim_File ' ' opts ' ' append_unix]); diff --git a/openEMS/matlab/Tutorials/CylindricalWave_CC.m b/openEMS/matlab/Tutorials/CylindricalWave_CC.m index d55c470..b81d2d4 100644 --- a/openEMS/matlab/Tutorials/CylindricalWave_CC.m +++ b/openEMS/matlab/Tutorials/CylindricalWave_CC.m @@ -28,7 +28,7 @@ exite_offset = 1300; excite_angle = 45; %% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%% -FDTD = InitFDTD(100000,1e-4,'CoordSystem',1,'MultiGrid',split); +FDTD = InitFDTD('NrTS', 100000, 'EndCriteria', 1e-4, 'CoordSystem', 1, 'MultiGrid', split); FDTD = SetGaussExcite(FDTD,f0,f0/2); BC = [0 3 0 0 0 0]; % pml in positive r-direction FDTD = SetBoundaryCond(FDTD,BC); diff --git a/openEMS/matlab/Tutorials/RadarUWBTutorial.m b/openEMS/matlab/Tutorials/RadarUWBTutorial.m new file mode 100644 index 0000000..6fd2f5b --- /dev/null +++ b/openEMS/matlab/Tutorials/RadarUWBTutorial.m @@ -0,0 +1,236 @@ +% Tutorial on time delay and signal integrity for radar +% and UWB applications +% +% Tested with +% - Octave 4.0 +% - openEMS v0.0.35 +% +% Author: Georg Michel, 2016 + + clear; + close all; + +physical_constants; + +% --- start of configuration section --- + +% In radar and ultrawideband applications it is important to know the +% delay and fidelity of RF pulses. The delay is the retardation of the +% signal from the source to the phase center of the antenna. It is +% composed out of linear delay, dispersion and minimum-phase +% delay. Dispersion due to waveguides or frequency-dependent +% permittivity and minimum-phase delay due to resonances will degrade +% the fidelity which is the normalized similarity between excitation and +% radiated signal. In this tutorial you can examine the performance of a +% simple ultrawideband (UWB) monopole. The delay and fidelity of this +% antenna are calculated and plotted. You can compare these properties +% in different channels. +% +% The Gaussian excitation is set to the same 3dB bandwidth as the +% channels of the IEEE 802.15.4 UWB PHY. One exeption is channel4twice +% which has the double bandwidth of channel 4. It can be seen that the +% delay is larger and the fidelity is smaller in the vicinity of the +% (undesired) resonances of the antenna. Note that for a real UWB system +% the total delay and fidelity result from both the transmitting and +% receiving antenna or twice the delay and the square of the fidelity +% for monostatic radar. +% +% The resolution of the delay will depend on the 'Oversampling' +% parameter to InitFDTD. See the description of DelayFidelity +% +% In the configuration section below you can uncomment the respective +% parameter settings. As an exercise, you can examine how the permittivity +% of the substrate influences gain, delay and fidelity. + + +%suffix = "channel1"; +%f_0 = 3.5e9; % center frequency of the channel +%f_c = 0.25e9 / 0.3925; % 3dB bandwidth is 0.3925 times 20dB bandwidth for Gaussian excitation + +%suffix = "channel2"; +%f_0 = 4.0e9; % center frequency of the channel +%f_c = 0.25e9 / 0.3925; + +%suffix = "channel3"; +%f_0 = 4.5e9; % center frequency of the channel +%f_c = 0.25e9 / 0.3925; + +suffix = "channel4"; +f_0 = 4.0e9; % center frequency of the channel +f_c = 0.5e9 / 0.3925; + +%suffix = "channel5"; +%f_0 = 6.5e9; % center frequency of the channel +%f_c = 0.25e9 / 0.3925; + +%suffix = "channel7"; +%f_0 = 6.5e9; % center frequency of the channel +%f_c = 0.5e9 / 0.3925; + +%suffix = "channel4twice"; % this is just to demonstrate the degradation of the fidelity with increasing bandwidth +%f_0 = 4.0e9; % center frequency of the channel +%f_c = 1e9 / 0.3925; + +tilt = 45 * pi / 180; % polarization tilt angle against co-polarization (90DEG is cross polarized) + +% --- end of configuration section --- + +% path and filename setup +Sim_Path = 'tmp'; +Sim_CSX = 'uwb.xml'; + +% properties of the substrate +substrate.epsR = 4; % FR4 +substrate.height = 0.707; +substrate.cells = 3; % thickness in cells + +% size of the monopole and the gap to the ground plane +gap = 0.62; % 0.5 +patchsize = 14; + +% we will use millimeters +unit = 1e-3; + +% set the resolution for the finer structures, e.g. the antenna gap +fineResolution = C0 / (f_0 + f_c) / sqrt(substrate.epsR) / unit / 40; +% set the resolution for the coarser structures, e.g. the surrounding air +coarseResolution = C0/(f_0 + f_c) / unit / 20; + + +% initialize the CSX structure +CSX = InitCSX(); + +% add the properties which are used to model the antenna +CSX = AddMetal(CSX, 'Ground' ); +CSX = AddMetal(CSX, 'Patch'); +CSX = AddMetal(CSX, 'Line'); +CSX = AddMaterial(CSX, 'Substrate' ); +CSX = SetMaterialProperty(CSX, 'Substrate', 'Epsilon', substrate.epsR); + +% define the supstrate and sheet-like primitives for the properties +CSX = AddBox(CSX, 'Substrate', 1, [-16, -16, -substrate.height], [16, 18, 0]); +CSX = AddBox(CSX, 'Ground', 2, [-16, -16, -substrate.height], [16, 0, -substrate.height]); +CSX = AddBox(CSX, 'Line', 2, [-1.15, -16, 0], [1.15, gap, 0]); +CSX = AddBox(CSX, 'Patch', 2, [-patchsize/2, gap, 0], [patchsize/2, gap + patchsize, 0]); + +% setup a mesh +mesh.x = []; +mesh.y = []; + +% two mesh lines for the metal coatings of teh substrate +mesh.z = linspace(-substrate.height, 0, substrate.cells +1); + +% find optimal mesh lines for the patch and ground, not yes the microstrip line +mesh = DetectEdges(CSX, mesh, 'SetProperty',{'Patch', 'Ground'}, '2D_Metal_Edge_Res', fineResolution/2); + +%replace gap mesh lines which are too close by a single mesh line +tooclose = find (diff(mesh.y) < fineResolution/4); +if ~isempty(tooclose) + mesh.y(tooclose) = (mesh.y(tooclose) + mesh.y(tooclose+1))/2; + mesh.y(tooclose + 1) = []; +endif + +% store the microstrip edges in a temporary variable +meshline = DetectEdges(CSX, [], 'SetProperty', 'Line', '2D_Metal_Edge_Res', fineResolution/2); +% as well as the edges of the substrate (without 1/3 - 2/3 rule) +meshsubstrate = DetectEdges(CSX, [], 'SetProperty', 'Substrate'); +% add only the x mesh lines of the microstrip +mesh.x = [mesh.x meshline.x]; +% and only the top of the substrate, the other edges are covered by the ground plane +mesh.y = [mesh.y, meshsubstrate.y(end)]; % top of substrate + +% for now we have only the edges, now calculate mesh lines inbetween +mesh = SmoothMesh(mesh, fineResolution); + +% add the outer boundary +mesh.x = [mesh.x -60, 60]; +mesh.y = [mesh.y, -60, 65]; +mesh.z = [mesh.z, -46, 45]; + +% add coarse mesh lines for the free space +mesh = SmoothMesh(mesh, coarseResolution); + +% define the grid +CSX = DefineRectGrid( CSX, unit, mesh); +% and the feeding port +[CSX, port] = AddLumpedPort( CSX, 999, 1, 50, [-1.15, meshline.y(2), -substrate.height], [1.15, meshline.y(2), 0], [0 0 1], true); + +%setup a NF2FF box for the calculation of the far field +start = [mesh.x(10) mesh.y(10) mesh.z(10)]; +stop = [mesh.x(end-9) mesh.y(end-9) mesh.z(end-9)]; +[CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop); + +% initialize the FDTD structure with excitation and open boundary conditions +FDTD = InitFDTD( 'NrTs', 30000, 'EndCriteria', 1e-5, 'OverSampling', 20); +FDTD = SetGaussExcite(FDTD, f_0, f_c ); +BC = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8'}; +FDTD = SetBoundaryCond(FDTD, BC ); + + +% remove old data, show structure, calculate new data +[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory +[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder + +% write the data to the working directory +WriteOpenEMS([Sim_Path '/' Sim_CSX], FDTD, CSX); +% show the geometry for checking +CSXGeomPlot([Sim_Path '/' Sim_CSX]); +% run the simulation +RunOpenEMS( Sim_Path, Sim_CSX); + +% plot amplitude and phase of the reflection coefficient +freq = linspace(f_0-f_c, f_0+f_c, 200); +port = calcPort(port, Sim_Path, freq); +s11 = port.uf.ref ./ port.uf.inc; +s11phase = unwrap(arg(s11)); +figure %("visible", "off"); % use this to plot only into files at the end of this script +ax = plotyy( freq/1e6, 20*log10(abs(s11)), freq/1e6, s11phase); +grid on +title( ['reflection coefficient ', suffix, ' S_{11}']); +xlabel( 'frequency f / MHz' ); +ylabel( ax(1), 'reflection coefficient |S_{11}|' ); +ylabel(ax(2), 'S_{11} phase (rad)'); + +% define an azimuthal trace around the monopole +phi = [0] * pi / 180; +theta = [-180:10:180] * pi / 180; + +% calculate the delay, the fidelity and the farfield +[delay, fidelity, nf2ff] = DelayFidelity(nf2ff, port, Sim_Path, sin(tilt), cos(tilt), theta, phi, f_0, f_c, 'Mode', 1); + +%plot the gain at (close to) f_0 +f_0_nearest_ind = nthargout(2, @min, abs(nf2ff.freq -f_0)); +%turn directivity into gain +nf2ff.Dmax(f_0_nearest_ind) *= nf2ff.Prad(f_0_nearest_ind) / calcPort(port, Sim_Path, nf2ff.freq(f_0_nearest_ind)).P_inc; +figure %("visible", "off"); +polarFF(nf2ff, 'xaxis', 'theta', 'freq_index', f_0_nearest_ind, 'logscale', [-4, 4]); +title(["gain ", suffix, " / dBi"]); + + +% We trick polarFF into plotting the delay in mm because +% the axes of the vanilla polar plot can not be scaled. +plotvar = delay * C0 * 1000; +maxplot = 80; +minplot = 30; +nf2ff.Dmax(1) = 10^(max(plotvar)/10); +nf2ff.E_norm{1} = 10.^(plotvar/20); +figure %("visible", "off"); +polarFF(nf2ff, 'xaxis', 'theta', 'logscale', [minplot, maxplot]); +title(["delay ", suffix, " / mm"]); + +% The same for the fidelity in percent. +plotvar = fidelity * 100; +maxplot = 100; +minplot = 98; +nf2ff.Dmax(1) = 10^(max(plotvar)/10); +nf2ff.E_norm{1} = 10.^(plotvar/20); +figure %("visible", "off"); +polarFF(nf2ff, 'xaxis', 'theta', 'logscale', [minplot, maxplot]); +title(["fidelity ", suffix, " / %"]); + +% save the plots in order to compare them afer simulating the different channels +print(1, ["s11_", suffix, ".png"]); +print(2, ["farfield_", suffix, ".png"]); +print(3, ["delay_mm_", suffix, ".png"]); +print(4, ["fidelity_", suffix, ".png"]); +return; \ No newline at end of file diff --git a/openEMS/matlab/Tutorials/Simple_Patch_Antenna.m b/openEMS/matlab/Tutorials/Simple_Patch_Antenna.m index ee2a8f0..1e24c4c 100644 --- a/openEMS/matlab/Tutorials/Simple_Patch_Antenna.m +++ b/openEMS/matlab/Tutorials/Simple_Patch_Antenna.m @@ -1,20 +1,20 @@ -% -% Tutorials / simple patch antenna +%% Simple Patch Antenna Tutorial % % Describtion at: -% http://openems.de/index.php/Tutorial:_Simple_Patch_Antenna +% % % Tested with % - Matlab 2013a / Octave 4.0 -% - openEMS v0.0.33 +% - openEMS v0.0.35 % -% (C) 2010-2015 Thorsten Liebig +% (C) 2010-2017 Thorsten Liebig +%% close all clear clc -%% setup the simulation +%% Setup the Simulation physical_constants; unit = 1e-3; % all length in mm @@ -38,7 +38,7 @@ feed.R = 50; %feed resistance % size of the simulation box SimBox = [200 200 150]; -%% setup FDTD parameter & excitation function +%% Setup FDTD Parameter & Excitation Function f0 = 2e9; % center frequency fc = 1e9; % 20 dB corner frequency FDTD = InitFDTD( 'NrTs', 30000 ); @@ -46,7 +46,7 @@ FDTD = SetGaussExcite( FDTD, f0, fc ); BC = {'MUR' 'MUR' 'MUR' 'MUR' 'MUR' 'MUR'}; % boundary conditions FDTD = SetBoundaryCond( FDTD, BC ); -%% setup CSXCAD geometry & mesh +%% Setup CSXCAD Geometry & Mesh CSX = InitCSX(); %initialize the mesh with the "air-box" dimensions @@ -54,13 +54,13 @@ mesh.x = [-SimBox(1)/2 SimBox(1)/2]; mesh.y = [-SimBox(2)/2 SimBox(2)/2]; mesh.z = [-SimBox(3)/3 SimBox(3)*2/3]; -%% create patch +% Create Patch CSX = AddMetal( CSX, 'patch' ); % create a perfect electric conductor (PEC) start = [-patch.width/2 -patch.length/2 substrate.thickness]; stop = [ patch.width/2 patch.length/2 substrate.thickness]; CSX = AddBox(CSX,'patch',10,start,stop); % add a box-primitive to the metal property 'patch' -%% create substrate +% Create Substrate CSX = AddMaterial( CSX, 'substrate' ); CSX = SetMaterialProperty( CSX, 'substrate', 'Epsilon', substrate.epsR, 'Kappa', substrate.kappa ); start = [-substrate.width/2 -substrate.length/2 0]; @@ -70,18 +70,18 @@ CSX = AddBox( CSX, 'substrate', 0, start, stop ); % add extra cells to discretize the substrate thickness mesh.z = [linspace(0,substrate.thickness,substrate.cells+1) mesh.z]; -%% create ground (same size as substrate) +% Create Ground same size as substrate CSX = AddMetal( CSX, 'gnd' ); % create a perfect electric conductor (PEC) start(3)=0; stop(3) =0; CSX = AddBox(CSX,'gnd',10,start,stop); -%% apply the excitation & resist as a current source +% Apply the Excitation & Resist as a Current Source start = [feed.pos 0 0]; stop = [feed.pos 0 substrate.thickness]; [CSX port] = AddLumpedPort(CSX, 5 ,1 ,feed.R, start, stop, [0 0 1], true); -%% finalize the mesh +% Finalize the Mesh % detect all edges except of the patch mesh = DetectEdges(CSX, mesh,'ExcludeProperty','patch'); % detect and set a special 2D metal edge mesh for the patch @@ -90,35 +90,41 @@ mesh = DetectEdges(CSX, mesh,'SetProperty','patch','2D_Metal_Edge_Res', c0/(f0+f mesh = SmoothMesh(mesh, c0/(f0+fc)/unit/20); CSX = DefineRectGrid(CSX, unit, mesh); -%% add a nf2ff calc box; size is 3 cells away from MUR boundary condition +CSX = AddDump(CSX,'Hf', 'DumpType', 11, 'Frequency',[2.4e9]); +CSX = AddBox(CSX,'Hf',10,[-substrate.width -substrate.length -10*substrate.thickness],[substrate.width +substrate.length 10*substrate.thickness]); %assign box + +% add a nf2ff calc box; size is 3 cells away from MUR boundary condition start = [mesh.x(4) mesh.y(4) mesh.z(4)]; stop = [mesh.x(end-3) mesh.y(end-3) mesh.z(end-3)]; [CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop); -%% prepare simulation folder +%% Prepare and Run Simulation Sim_Path = 'tmp_Patch_Ant'; Sim_CSX = 'patch_ant.xml'; +% create an empty working directory [status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory [status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder -%% write openEMS compatible xml-file +% write openEMS compatible xml-file WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX ); -%% show the structure +% show the structure CSXGeomPlot( [Sim_Path '/' Sim_CSX] ); -%% run openEMS +% run openEMS RunOpenEMS( Sim_Path, Sim_CSX); -%% postprocessing & do the plots +%% Postprocessing & Plots freq = linspace( max([1e9,f0-fc]), f0+fc, 501 ); port = calcPort(port, Sim_Path, freq); -Zin = port.uf.tot ./ port.if.tot; -s11 = port.uf.ref ./ port.uf.inc; +%% Smith chart port reflection +plotRefl(port, 'threshold', -10) +title( 'reflection coefficient' ); % plot feed point impedance +Zin = port.uf.tot ./ port.if.tot; figure plot( freq/1e6, real(Zin), 'k-', 'Linewidth', 2 ); hold on @@ -130,6 +136,7 @@ ylabel( 'impedance Z_{in} / Ohm' ); legend( 'real', 'imag' ); % plot reflection coefficient S11 +s11 = port.uf.ref ./ port.uf.inc; figure plot( freq/1e6, 20*log10(abs(s11)), 'k-', 'Linewidth', 2 ); grid on @@ -139,7 +146,7 @@ ylabel( 'reflection coefficient |S_{11}|' ); drawnow -%% NFFF contour plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% NFFF Plots %find resonance frequncy from s11 f_res_ind = find(s11==min(s11)); f_res = freq(f_res_ind); @@ -166,7 +173,7 @@ plotFFdB(nf2ff,'xaxis','theta','param',[1 2]) drawnow -%% +% Show 3D pattern disp( 'calculating 3D far field pattern and dumping to vtk (use Paraview to visualize)...' ); thetaRange = (0:2:180); phiRange = (0:2:360) - 180; diff --git a/openEMS/matlab/Tutorials/StripLine2MSL.m b/openEMS/matlab/Tutorials/StripLine2MSL.m new file mode 100644 index 0000000..64eeaa1 --- /dev/null +++ b/openEMS/matlab/Tutorials/StripLine2MSL.m @@ -0,0 +1,133 @@ +% +% Stripline to Microstrip Line Transition +% +% Describtion at: +% +% +% Tested with +% - Octave 4.0 +% - openEMS v0.0.35 +% +% (C) 2017 Thorsten Liebig + +close all +clear +clc + +%% Setup the Simulation +physical_constants; +unit = 1e-6; % specify everything in um + +line_length = 15000; % line length of strip line and microstrip line +substrate_width = 6000; +air_spacer = 4000; % air spacer above the substrate + +msl_width = 500; +msl_substrate_thickness = 254; + +strip_width = 500; +strip_substrate_thickness = 512; + +connect_via_rad = 500/2; +connect_via_gap = 1250/2; + +substrate_epr = 3.66; +substrate_kappa = 1e-3 * 2*pi*2.45e9 * EPS0*substrate_epr; % substrate losses + +f_max = 10e9; +resolution = 250; +edge_res = 25; +feed_shift = 2500; +meas_shift = 5000; + +%% Setup FDTD Parameters & Excitation Function +FDTD = InitFDTD(); +FDTD = SetGaussExcite( FDTD, f_max/2, f_max/2); +BC = {'PML_8' 'PML_8' 'PEC' 'PEC' 'PEC' 'MUR'}; +FDTD = SetBoundaryCond( FDTD, BC ); + +%% Setup CSXCAD Geometry & Mesh +CSX = InitCSX(); +edge_mesh = [-1/3 2/3]*edge_res; % 1/3 - 2/3 rule for 2D metal edges + +mesh.x = SmoothMeshLines( [-connect_via_gap 0 connect_via_gap], 2*edge_res, 1.5 ); +mesh.x = SmoothMeshLines( [-line_length mesh.x line_length], resolution, 1.5); +mesh.y = SmoothMeshLines( [0 msl_width/2+edge_mesh substrate_width/2], resolution/4 , 1.5); +mesh.y = sort(unique([-mesh.y mesh.y])); +mesh.z = SmoothMeshLines( [linspace(-strip_substrate_thickness,0,5) linspace(0,strip_substrate_thickness,5) linspace(strip_substrate_thickness,msl_substrate_thickness+strip_substrate_thickness,5) 2*strip_substrate_thickness+air_spacer] , resolution ); +CSX = DefineRectGrid( CSX, unit, mesh ); + +% Create Substrate +CSX = AddMaterial( CSX, 'RO4350B' ); +CSX = SetMaterialProperty( CSX, 'RO4350B', 'Epsilon', substrate_epr, 'Kappa', substrate_kappa ); +start = [mesh.x(1), mesh.y(1), -strip_substrate_thickness]; +stop = [mesh.x(end), mesh.y(end), +strip_substrate_thickness+msl_substrate_thickness]; +CSX = AddBox( CSX, 'RO4350B', 0, start, stop ); + +% Create a PEC called 'metal' and 'gnd' +CSX = AddMetal( CSX, 'gnd' ); +CSX = AddMetal( CSX, 'metal' ); + +% Create strip line port (incl. metal stip line) +start = [-line_length -strip_width/2 0]; +stop = [0 +strip_width/2 0]; +[CSX,port{1}] = AddStripLinePort( CSX, 100, 1, 'metal', start, stop, strip_substrate_thickness, 'x', [0 0 -1], 'ExcitePort', true, 'FeedShift', feed_shift, 'MeasPlaneShift', meas_shift ); + +% Create MSL port on top +start = [line_length -strip_width/2 strip_substrate_thickness+msl_substrate_thickness]; +stop = [0 +strip_width/2 strip_substrate_thickness]; +[CSX,port{2}] = AddMSLPort( CSX, 100, 2, 'metal', start, stop, 'x', [0 0 -1], 'MeasPlaneShift', meas_shift ); + +% transitional via +start = [0, 0, 0]; +stop = [0, 0, strip_substrate_thickness+msl_substrate_thickness]; +CSX = AddCylinder(CSX, 'metal', 100, start, stop, connect_via_rad); + +% metal plane between strip line and MSL, including hole for transition +p(1,1) = mesh.x(1); +p(2,1) = mesh.y(1); +p(1,2) = 0; +p(2,2) = mesh.y(1); +for a = linspace(-pi, pi, 21) + p(1,end+1) = connect_via_gap*sin(a); + p(2,end) = connect_via_gap*cos(a); +endfor +p(1,end+1) = 0; +p(2,end ) = mesh.y(1); +p(1,end+1) = mesh.x(end); +p(2,end ) = mesh.y(1); +p(1,end+1) = mesh.x(end); +p(2,end ) = mesh.y(end); +p(1,end+1) = mesh.x(1); +p(2,end ) = mesh.y(end); +CSX = AddPolygon( CSX, 'gnd', 1, 'z', strip_substrate_thickness, p); + +%% Write/Show/Run the openEMS compatible xml-file +Sim_Path = 'tmp'; +Sim_CSX = 'strip2msl.xml'; + +[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory +[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder + +WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX ); +CSXGeomPlot( [Sim_Path '/' Sim_CSX] ); +RunOpenEMS( Sim_Path, Sim_CSX ); + +%% Post-Processing +close all +f = linspace( 0, f_max, 1601 ); +port = calcPort( port, Sim_Path, f, 'RefImpedance', 50); + +s11 = port{1}.uf.ref./ port{1}.uf.inc; +s21 = port{2}.uf.ref./ port{1}.uf.inc; + +plot(f/1e9,20*log10(abs(s11)),'k-','LineWidth',2); +hold on; +grid on; +plot(f/1e9,20*log10(abs(s21)),'r--','LineWidth',2); +legend('S_{11}','S_{21}'); +ylabel('S-Parameter (dB)','FontSize',12); +xlabel('frequency (GHz) \rightarrow','FontSize',12); +ylim([-40 2]); + + diff --git a/openEMS/matlab/h5readatt_octave.cc b/openEMS/matlab/h5readatt_octave.cc index 8bd58d0..13e1765 100755 --- a/openEMS/matlab/h5readatt_octave.cc +++ b/openEMS/matlab/h5readatt_octave.cc @@ -5,13 +5,7 @@ // this special treatment is necessary because Win32-Octave ships with a very old hdf5 version (1.6.10) void CloseH5Object(hid_t obj) { -#if ((H5_VERS_MAJOR == 1) && (H5_VERS_MINOR == 6)) - // try group close, than Dataset close - if (H5Gclose(obj)<0) - H5Dclose(obj); -#else H5Oclose(obj); -#endif } DEFUN_DLD (h5readatt_octave, args, nargout, "h5readatt_octave(,,)") @@ -30,7 +24,7 @@ DEFUN_DLD (h5readatt_octave, args, nargout, "h5readatt_octave(,,, threshold; +lowerind = s11dB(1:end-1) > threshold & s11dB(2:end) < threshold; +minind = nthargout(2, @min, s11dB); +handle1 = plot(s11(lowerind),['<','b']); +handle2 = plot(s11(upperind),['>','b']); +handle3 = plot(s11(minind),['*', 'b']); +llegend = num2str(port.f(lowerind)(1)/1e6, ffmt); +ulegend = num2str(port.f(upperind)(1)/1e6, ffmt); + +if nnz(lowerind) > 1 + for i= 2:nnz(lowerind) + llegend = strjoin({llegend, num2str(port.f(lowerind)(i)/1e6, ffmt)}, ', '); + end +end + +if nnz(upperind) > 1 + for i= 2:nnz(upperind) + ulegend = strjoin({ulegend, num2str(port.f(upperind)(i)/1e6, ffmt)}, ', '); + end +end + +legend([handle1, handle2, handle3], {[llegend, " MHz"], ... + [ulegend, " MHz"], ... + [num2str(20*log10(abs(s11(minind))), "%4.0f"), ... + "dB @ ", num2str(port.f(minind)/1e6, ffmt), " MHz"]}); +h = plot(s11); + +if (nargout == 0) + clear h; +end + +end \ No newline at end of file diff --git a/openEMS/matlab/setup.m b/openEMS/matlab/setup.m index df56b25..e208ad4 100644 --- a/openEMS/matlab/setup.m +++ b/openEMS/matlab/setup.m @@ -5,7 +5,7 @@ function setup() % % openEMS matlab/octave interface % ----------------------- -% author: Thorsten Liebig (2011) +% author: Thorsten Liebig (2011-2017) disp('setting up openEMS matlab/octave interface') @@ -16,18 +16,22 @@ cd(dir); if isOctave() disp('compiling oct files') - fflush(stdout) + fflush(stdout); if isunix - [res, fn] = unix('find /usr/lib -name libhdf5.so'); - if length(fn)>0 - [hdf5lib_dir, hdf5lib_fn] = fileparts(fn); + [res, fn_so] = unix('find /usr/lib -name libhdf5.so'); + [res, fn_h] = unix('find /usr/include -name hdf5.h'); + if length(fn_so)>0 && length(fn_h)>0 + [hdf5lib_dir, hdf5lib_fn] = fileparts(fn_so); disp(["HDF5 library path found at: " hdf5lib_dir]) - mkoctfile(["-L" hdf5lib_dir ],"-lhdf5 -DH5_USE_16_API", "h5readatt_octave.cc") + + [hdf5inc_dir, hdf5inc_fn] = fileparts(fn_h); + disp(["HDF5 include path found at: " hdf5inc_dir]) + mkoctfile(["-L" hdf5lib_dir " -I" hdf5inc_dir],"-lhdf5", "h5readatt_octave.cc") else - mkoctfile -lhdf5 -DH5_USE_16_API h5readatt_octave.cc + mkoctfile -lhdf5 h5readatt_octave.cc end else - mkoctfile -lhdf5 -DH5_USE_16_API h5readatt_octave.cc + mkoctfile -lhdf5 h5readatt_octave.cc end else disp('Matlab does not need this function. It is Octave only.') diff --git a/openEMS/nf2ff/CMakeLists.txt b/openEMS/nf2ff/CMakeLists.txt index d43d519..86811b4 100644 --- a/openEMS/nf2ff/CMakeLists.txt +++ b/openEMS/nf2ff/CMakeLists.txt @@ -38,6 +38,7 @@ TARGET_LINK_LIBRARIES( nf2ff tinyxml ${HDF5_LIBRARIES} ${Boost_LIBRARIES} + ${MPI_LIBRARIES} ) ADD_EXECUTABLE( nf2ff_bin main.cpp ) @@ -47,5 +48,7 @@ TARGET_LINK_LIBRARIES(nf2ff_bin nf2ff) INSTALL(TARGETS nf2ff_bin DESTINATION bin) INSTALL(TARGETS nf2ff DESTINATION lib${LIB_SUFFIX}) +INSTALL(FILES ${HEADERS} DESTINATION include/openEMS) + #TODO tarball, debug, release diff --git a/openEMS/openEMS_MPI.pro b/openEMS/openEMS_MPI.pro deleted file mode 100644 index d564c2b..0000000 --- a/openEMS/openEMS_MPI.pro +++ /dev/null @@ -1,5 +0,0 @@ -# enable MPI support -!win32:CONFIG += MPI_SUPPORT - -include(openEMS.pro) - diff --git a/openEMS/openems.cpp b/openEMS/openems.cpp index 72476d3..1d08ef7 100644 --- a/openEMS/openems.cpp +++ b/openEMS/openems.cpp @@ -418,11 +418,13 @@ bool openEMS::SetupProcessing() } if (CylinderCoords) proc->SetMeshType(Processing::CYLINDRICAL_MESH); - if ((pb->GetProbeType()==1) || (pb->GetProbeType()==3) || (pb->GetProbeType()==11)) + if ((pb->GetProbeType()==1) || (pb->GetProbeType()==3)) { proc->SetDualTime(true); proc->SetDualMesh(true); } + if (pb->GetProbeType()==11) + proc->SetDualTime(true); proc->SetProcessInterval(Nyquist/m_OverSampling); if (pb->GetStartTime()>0 || pb->GetStopTime()>0) proc->SetProcessStartStopTime(pb->GetStartTime(), pb->GetStopTime()); @@ -715,9 +717,7 @@ bool openEMS::Parse_XML_FDTDSetup(TiXmlElement* FDTD_Opts) ihelp = 0; FDTD_Opts->QueryIntAttribute("OverSampling",&ihelp); - if (ihelp<2) - this->SetOverSampling(2); - else + if (ihelp>1) this->SetOverSampling(ihelp); // check for cell constant material averaging @@ -821,6 +821,7 @@ bool openEMS::Parse_XML_FDTDSetup(TiXmlElement* FDTD_Opts) this->SetTimeStep(dhelp); if (FDTD_Opts->QueryDoubleAttribute("TimeStepFactor",&dhelp)==TIXML_SUCCESS) this->SetTimeStepFactor(dhelp); + return true; } void openEMS::SetGaussExcite(double f0, double fc) diff --git a/openEMS/tools/hdf5_file_reader.cpp b/openEMS/tools/hdf5_file_reader.cpp index 5e81614..a6e1359 100644 --- a/openEMS/tools/hdf5_file_reader.cpp +++ b/openEMS/tools/hdf5_file_reader.cpp @@ -68,7 +68,7 @@ bool HDF5_File_Reader::OpenGroup(hid_t &file, hid_t &group, string groupName) return 0; } - group = H5Gopen(file, groupName.c_str() ); + group = H5Gopen(file, groupName.c_str(), H5P_DEFAULT ); if (group<0) { cerr << "HDF5_File_Reader::OpenGroup: can't open group """ << groupName << """" << endl; @@ -198,7 +198,7 @@ bool HDF5_File_Reader::ReadDataSet(string ds_name, hsize_t &nDim, hsize_t* &dims return false; } - hid_t dataset = H5Dopen(hdf5_file, ds_name.c_str() ); + hid_t dataset = H5Dopen(hdf5_file, ds_name.c_str(), H5P_DEFAULT ); if (dataset<0) { cerr << "HDF5_File_Reader::ReadDataSet: dataset not found" << endl; diff --git a/openEMS/tools/hdf5_file_writer.cpp b/openEMS/tools/hdf5_file_writer.cpp index 81c5692..e3a2113 100644 --- a/openEMS/tools/hdf5_file_writer.cpp +++ b/openEMS/tools/hdf5_file_writer.cpp @@ -52,7 +52,7 @@ hid_t HDF5_File_Writer::OpenGroup(hid_t hdf5_file, string group) vector results; boost::split(results, group, boost::is_any_of("/")); - hid_t grp=H5Gopen(hdf5_file,"/"); + hid_t grp=H5Gopen(hdf5_file,"/", H5P_DEFAULT); if (grp<0) { cerr << "HDF5_File_Writer::OpenGroup: Error, opening root group " << endl; @@ -66,7 +66,7 @@ hid_t HDF5_File_Writer::OpenGroup(hid_t hdf5_file, string group) { if (H5Lexists(grp, results.at(n).c_str(), H5P_DEFAULT)) { - grp = H5Gopen(grp, results.at(n).c_str()); + grp = H5Gopen(grp, results.at(n).c_str(), H5P_DEFAULT); H5Gclose(old_grp); if (grp<0) { @@ -76,7 +76,7 @@ hid_t HDF5_File_Writer::OpenGroup(hid_t hdf5_file, string group) } else { - grp = H5Gcreate(grp,results.at(n).c_str(),0); + grp = H5Gcreate(grp,results.at(n).c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); H5Gclose(old_grp); if (grp<0) { @@ -140,7 +140,7 @@ bool HDF5_File_Writer::WriteRectMesh(unsigned int const* numLines, float const* return false; } - hid_t mesh_grp = H5Gcreate(hdf5_file,"/Mesh",0); + hid_t mesh_grp = H5Gcreate(hdf5_file,"/Mesh", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); if (mesh_grp<0) { cerr << "HDF5_File_Writer::WriteRectMesh: Error, creating group ""/Mesh"" failed" << endl; @@ -165,7 +165,7 @@ bool HDF5_File_Writer::WriteRectMesh(unsigned int const* numLines, float const* { hsize_t dims[1]={numLines[n]}; hid_t space = H5Screate_simple(1, dims, NULL); - hid_t dataset = H5Dcreate(mesh_grp, names[n].c_str(), H5T_NATIVE_FLOAT, space, H5P_DEFAULT); + hid_t dataset = H5Dcreate(mesh_grp, names[n].c_str(), H5T_NATIVE_FLOAT, space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); float* array = new float[numLines[n]]; for (unsigned int i=0; i