From 3003a81929cc3eef567fc77f7ada75bf8da5413d Mon Sep 17 00:00:00 2001 From: Jan Behrens Date: Fri, 29 May 2020 16:45:36 -0400 Subject: [PATCH] Update to Kassiopeia 3.7.3 --- KEMField/CMakeLists.txt | 2 +- KEMField/Source/Bindings/CMakeLists.txt | 2 + .../src/KElectricZHFieldSolverBuilder.cc | 2 +- ...HarmonicMagnetostaticFieldSolverBuilder.cc | 2 +- .../KElectrostaticLinearFieldBuilder.hh | 55 +++ .../KElectrostaticPotentialmapBuilder.hh | 4 + .../src/KElectrostaticLinearFieldBuilder.cc | 22 + .../src/KElectrostaticPotentialmapBuilder.cc | 1 + .../include/KMagnetostaticFieldmapBuilder.hh | 4 + .../include/KRampedMagneticFieldBuilder.hh | 3 + .../src/KMagnetostaticFieldmapBuilder.cc | 1 + .../include/KZonalHarmonicContainer.hh | 43 ++ .../src/KZHCoefficientGeneratorElement.cc | 11 +- .../KElectromagnetZonalHarmonicFieldSolver.hh | 3 + .../KElectrostaticZonalHarmonicFieldSolver.hh | 3 + .../KElectromagnetZonalHarmonicFieldSolver.cc | 40 +- .../KElectrostaticZonalHarmonicFieldSolver.cc | 26 + .../include/KElectricZHFieldSolver.hh | 5 + .../Electric/src/KElectricZHFieldSolver.cc | 14 +- .../KZonalHarmonicMagnetostaticFieldSolver.hh | 9 + .../KZonalHarmonicMagnetostaticFieldSolver.cc | 21 + .../Interface/Fields/Electric/CMakeLists.txt | 2 + .../include/KElectrostaticLinearField.hh | 50 ++ .../Electric/src/KElectrostaticLinearField.cc | 89 ++++ .../KGeoBag/include/KGBEMConverter.hh | 6 +- .../Interface/KGeoBag/src/KGBEMConverter.cc | 28 +- .../KGeoBag/src/KGElectromagnetConverter.cc | 21 + .../src/KGElectrostaticBoundaryField.cc | 7 +- .../KGeoBag/src/KGStaticElectromagnetField.cc | 6 +- .../include/KElectrostaticPotentialmap.hh | 5 + .../include/KMagnetostaticFieldmap.hh | 6 + .../src/KElectrostaticPotentialmap.cc | 135 +++--- .../VTKPart2/src/KMagnetostaticFieldmap.cc | 107 ++++- KGeoBag/CMakeLists.txt | 2 +- KGeoBag/Source/Application/Source/KGeoBag.cc | 4 +- .../Include/KGROOTGeometryPainterBuilder.hh | 8 + .../Source/KGROOTGeometryPainterBuilder.cc | 2 + KGeoBag/Source/Core/Include/KGSpace.hh | 2 + KGeoBag/Source/Core/Include/KGSurface.hh | 1 + KGeoBag/Source/Core/Source/KGSpace.cc | 11 + KGeoBag/Source/Core/Source/KGSurface.cc | 10 + .../Root/Include/KGROOTGeometryPainter.hh | 8 +- .../Root/Source/KGROOTGeometryPainter.cc | 156 +++++- .../FieldCalculation/CMakeLists.txt | 2 + ...impleElectricFieldCalculatorAlongZaxis.cxx | 3 +- ...mpleElectricFieldCalculatorOverXYplane.cxx | 99 ++++ ...mpleMagneticFieldCalculatorOverXYplane.cxx | 94 ++++ .../Applications/Simulation/CMakeLists.txt | 19 +- .../Simulation/Source/Kassiopeia.cxx | 6 +- .../Simulation/Source/ParticleGenerator.cxx | 149 ++++++ Kassiopeia/Bindings/CMakeLists.txt | 2 + .../Simulation/Include/KSRootBuilder.h | 10 +- .../KSROOTZonalHarmonicsPainterBuilder.h | 109 +++++ .../KSROOTZonalHarmonicsPainterBuilder.cxx | 39 ++ Kassiopeia/CMakeLists.txt | 2 +- Kassiopeia/Visualization/CMakeLists.txt | 2 + .../Include/KSROOTZonalHarmonicsPainter.h | 88 ++++ .../Source/KSROOTZonalHarmonicsPainter.cxx | 453 ++++++++++++++++++ Kassiopeia/XML/CMakeLists.txt | 5 + .../XML/Examples/AnalyticSimulation.xml | 41 +- .../DipoleTrapMeshedSpaceSimulation.xml | 210 ++++---- .../XML/Examples/DipoleTrapSimulation.xml | 153 +++--- Kassiopeia/XML/Examples/Fields.xml | 326 +++++++++++++ Kassiopeia/XML/Examples/Generators.xml | 316 ++++++++++++ .../PhotoMultiplierTubeSimulation.xml | 351 ++++++++------ .../XML/Examples/QuadrupoleTrapSimulation.xml | 109 ++--- .../XML/Examples/ToricTrapSimulation.xml | 101 ++-- .../XML/Examples/VTK/AnalyticSimulation.xml | 115 ++--- .../VTK/DipoleTrapMeshedSpaceSimulation.xml | 264 +++++----- .../XML/Examples/VTK/DipoleTrapSimulation.xml | 233 ++++----- .../VTK/PhotoMultiplierTubeSimulation.xml | 409 +++++++++------- .../Examples/VTK/QuadrupoleTrapSimulation.xml | 161 +++---- .../XML/Examples/VTK/ToricTrapSimulation.xml | 147 +++--- Kassiopeia/XML/Validation/TestRampedField.xml | 57 +-- Kommon/Core/Initialization/KElementBase.cc | 17 +- Kommon/Core/Initialization/KXMLInitializer.cc | 9 +- Kommon/Root/Utility/KROOTWindow.cxx | 3 + 77 files changed, 3785 insertions(+), 1258 deletions(-) create mode 100644 KEMField/Source/Bindings/Fields/Electric/include/KElectrostaticLinearFieldBuilder.hh create mode 100644 KEMField/Source/Bindings/Fields/Electric/src/KElectrostaticLinearFieldBuilder.cc create mode 100644 KEMField/Source/Interface/Fields/Electric/include/KElectrostaticLinearField.hh create mode 100644 KEMField/Source/Interface/Fields/Electric/src/KElectrostaticLinearField.cc create mode 100644 Kassiopeia/Applications/FieldCalculation/Source/SimpleElectricFieldCalculatorOverXYplane.cxx create mode 100644 Kassiopeia/Applications/FieldCalculation/Source/SimpleMagneticFieldCalculatorOverXYplane.cxx create mode 100644 Kassiopeia/Applications/Simulation/Source/ParticleGenerator.cxx create mode 100644 Kassiopeia/Bindings/Visualization/Include/KSROOTZonalHarmonicsPainterBuilder.h create mode 100644 Kassiopeia/Bindings/Visualization/Source/KSROOTZonalHarmonicsPainterBuilder.cxx create mode 100644 Kassiopeia/Visualization/Include/KSROOTZonalHarmonicsPainter.h create mode 100644 Kassiopeia/Visualization/Source/KSROOTZonalHarmonicsPainter.cxx create mode 100644 Kassiopeia/XML/Examples/Fields.xml create mode 100644 Kassiopeia/XML/Examples/Generators.xml diff --git a/KEMField/CMakeLists.txt b/KEMField/CMakeLists.txt index 5d1754d0b..de47bc31c 100644 --- a/KEMField/CMakeLists.txt +++ b/KEMField/CMakeLists.txt @@ -7,7 +7,7 @@ endif() # KEMField version set(MODULE_VERSION_MAJOR 3) set(MODULE_VERSION_MINOR 7) -set(MODULE_VERSION_PATCH 0) +set(MODULE_VERSION_PATCH 3) set(MODULE_VERSION "${MODULE_VERSION_MAJOR}.${MODULE_VERSION_MINOR}.${MODULE_VERSION_PATCH}") project( KEMField ) diff --git a/KEMField/Source/Bindings/CMakeLists.txt b/KEMField/Source/Bindings/CMakeLists.txt index 90c3df8c8..83489fe69 100644 --- a/KEMField/Source/Bindings/CMakeLists.txt +++ b/KEMField/Source/Bindings/CMakeLists.txt @@ -22,6 +22,7 @@ add_cflag (KEMFIELD_USE_KOMMON_BINDINGS) ${CMAKE_CURRENT_SOURCE_DIR}/Fields/Electric/include/KBoundaryElementInfoDisplayBuilder.hh ${CMAKE_CURRENT_SOURCE_DIR}/Fields/Electric/include/KElectricQuadrupoleFieldBuilder.hh ${CMAKE_CURRENT_SOURCE_DIR}/Fields/Electric/include/KElectrostaticConstantFieldBuilder.hh + ${CMAKE_CURRENT_SOURCE_DIR}/Fields/Electric/include/KElectrostaticLinearFieldBuilder.hh ${CMAKE_CURRENT_SOURCE_DIR}/Fields/Electric/include/KElectrostaticBoundaryFieldBuilder.hh ${CMAKE_CURRENT_SOURCE_DIR}/Fields/Electric/include/KElectrostaticBoundaryFieldTimerBuilder.hh ${CMAKE_CURRENT_SOURCE_DIR}/Fields/Electric/include/KVTKViewerVisitorBuilder.hh @@ -64,6 +65,7 @@ add_cflag (KEMFIELD_USE_KOMMON_BINDINGS) ${CMAKE_CURRENT_SOURCE_DIR}/Fields/Electric/src/KBoundaryElementInfoDisplayBuilder.cc ${CMAKE_CURRENT_SOURCE_DIR}/Fields/Electric/src/KElectricQuadrupoleFieldBuilder.cc ${CMAKE_CURRENT_SOURCE_DIR}/Fields/Electric/src/KElectrostaticConstantFieldBuilder.cc + ${CMAKE_CURRENT_SOURCE_DIR}/Fields/Electric/src/KElectrostaticLinearFieldBuilder.cc ${CMAKE_CURRENT_SOURCE_DIR}/Fields/Electric/src/KElectrostaticBoundaryFieldBuilder.cc ${CMAKE_CURRENT_SOURCE_DIR}/Fields/Electric/src/KElectrostaticBoundaryFieldTimerBuilder.cc ${CMAKE_CURRENT_SOURCE_DIR}/Fields/Electric/src/KVTKViewerVisitorBuilder.cc diff --git a/KEMField/Source/Bindings/FieldSolvers/Electric/src/KElectricZHFieldSolverBuilder.cc b/KEMField/Source/Bindings/FieldSolvers/Electric/src/KElectricZHFieldSolverBuilder.cc index 0c5edc3e0..0ff3c94cf 100644 --- a/KEMField/Source/Bindings/FieldSolvers/Electric/src/KElectricZHFieldSolverBuilder.cc +++ b/KEMField/Source/Bindings/FieldSolvers/Electric/src/KElectricZHFieldSolverBuilder.cc @@ -17,7 +17,7 @@ namespace katrin template<> KElectricZHFieldSolverBuilder::~KComplexElement() {} STATICINT sKElectricZHFieldSolverBuilderStructure = - KElectricZHFieldSolverBuilder::Attribute("number_of_bifurcations") + + KElectricZHFieldSolverBuilder::Attribute("number_of_bifurcations") + KElectricZHFieldSolverBuilder::Attribute("convergence_ratio") + KElectricZHFieldSolverBuilder::Attribute("proximity_to_sourcepoint") + KElectricZHFieldSolverBuilder::Attribute("convergence_parameter") + diff --git a/KEMField/Source/Bindings/FieldSolvers/Magnetic/src/KZonalHarmonicMagnetostaticFieldSolverBuilder.cc b/KEMField/Source/Bindings/FieldSolvers/Magnetic/src/KZonalHarmonicMagnetostaticFieldSolverBuilder.cc index de03b3bc7..d79628f1b 100644 --- a/KEMField/Source/Bindings/FieldSolvers/Magnetic/src/KZonalHarmonicMagnetostaticFieldSolverBuilder.cc +++ b/KEMField/Source/Bindings/FieldSolvers/Magnetic/src/KZonalHarmonicMagnetostaticFieldSolverBuilder.cc @@ -17,7 +17,7 @@ namespace katrin template<> KZonalHarmonicMagnetostaticFieldSolverBuilder::~KComplexElement() {} STATICINT sKZonalHarmonicMagnetostaticFieldSolverStructure = - KZonalHarmonicMagnetostaticFieldSolverBuilder::Attribute("number_of_bifurcations") + + KZonalHarmonicMagnetostaticFieldSolverBuilder::Attribute("number_of_bifurcations") + KZonalHarmonicMagnetostaticFieldSolverBuilder::Attribute("convergence_ratio") + KZonalHarmonicMagnetostaticFieldSolverBuilder::Attribute("proximity_to_sourcepoint") + KZonalHarmonicMagnetostaticFieldSolverBuilder::Attribute("convergence_parameter") + diff --git a/KEMField/Source/Bindings/Fields/Electric/include/KElectrostaticLinearFieldBuilder.hh b/KEMField/Source/Bindings/Fields/Electric/include/KElectrostaticLinearFieldBuilder.hh new file mode 100644 index 000000000..71e763957 --- /dev/null +++ b/KEMField/Source/Bindings/Fields/Electric/include/KElectrostaticLinearFieldBuilder.hh @@ -0,0 +1,55 @@ +#ifndef KELECTROSTATICLINEARFIELDBUILDER_HH_ +#define KELECTROSTATICLINEARFIELDBUILDER_HH_ + +#include "KComplexElement.hh" +#include "KEMStreamableThreeVector.hh" +#include "KElectrostaticLinearField.hh" + +#include + +namespace katrin +{ + +typedef KComplexElement KElectrostaticLinearFieldBuilder; + +template<> inline bool KElectrostaticLinearFieldBuilder::AddAttribute(KContainer* aContainer) +{ + if (aContainer->GetName() == "name") { + std::string name; + aContainer->CopyTo(name); + this->SetName(name); + fObject->SetName(name); + return true; + } + if (aContainer->GetName() == "U1") { + aContainer->CopyTo(fObject, &KEMField::KElectrostaticLinearField::SetPotential1); + return true; + } + if (aContainer->GetName() == "U2") { + aContainer->CopyTo(fObject, &KEMField::KElectrostaticLinearField::SetPotential2); + return true; + } + if (aContainer->GetName() == "z1") { + aContainer->CopyTo(fObject, &KEMField::KElectrostaticLinearField::SetZ1); + return true; + } + if (aContainer->GetName() == "z2") { + aContainer->CopyTo(fObject, &KEMField::KElectrostaticLinearField::SetZ2); + return true; + } + if (aContainer->GetName() == "surface") { + // TODO + return false; + } + return false; +} + +template<> inline bool KElectrostaticLinearFieldBuilder::AddElement(KContainer* /*aContainer*/) +{ + return false; +} + +} // namespace katrin + + +#endif /* KELECTROSTATICCONSTANTFIELDBUILDER_HH_ */ diff --git a/KEMField/Source/Bindings/Fields/Electric/include/KElectrostaticPotentialmapBuilder.hh b/KEMField/Source/Bindings/Fields/Electric/include/KElectrostaticPotentialmapBuilder.hh index d846efc50..e23bb4b06 100644 --- a/KEMField/Source/Bindings/Fields/Electric/include/KElectrostaticPotentialmapBuilder.hh +++ b/KEMField/Source/Bindings/Fields/Electric/include/KElectrostaticPotentialmapBuilder.hh @@ -69,6 +69,10 @@ template<> inline bool KElectrostaticPotentialmapCalculatorBuilder::AddAttribute aContainer->CopyTo(fObject, &KEMField::KElectrostaticPotentialmapCalculator::SetForceUpdate); return true; } + if (aContainer->GetName() == "compute_field") { + aContainer->CopyTo(fObject, &KEMField::KElectrostaticPotentialmapCalculator::SetComputeField); + return true; + } if (aContainer->GetName() == "center") { KEMField::KEMStreamableThreeVector center; aContainer->CopyTo(center); diff --git a/KEMField/Source/Bindings/Fields/Electric/src/KElectrostaticLinearFieldBuilder.cc b/KEMField/Source/Bindings/Fields/Electric/src/KElectrostaticLinearFieldBuilder.cc new file mode 100644 index 000000000..e513ff849 --- /dev/null +++ b/KEMField/Source/Bindings/Fields/Electric/src/KElectrostaticLinearFieldBuilder.cc @@ -0,0 +1,22 @@ +#include "KElectrostaticLinearFieldBuilder.hh" + +#include "KEMToolboxBuilder.hh" + + +using namespace KEMField; +namespace katrin +{ + +template<> KElectrostaticLinearFieldBuilder::~KComplexElement() {} + +STATICINT sKEMToolBoxBuilder = KEMToolboxBuilder::ComplexElement("linear_electric_field"); + +STATICINT sKElectrostaticConstantFieldBuilder = + KElectrostaticLinearFieldBuilder::Attribute("name") + + KElectrostaticLinearFieldBuilder::Attribute("U1") + + KElectrostaticLinearFieldBuilder::Attribute("U2") + + KElectrostaticLinearFieldBuilder::Attribute("z1") + + KElectrostaticLinearFieldBuilder::Attribute("z2") + + KElectrostaticLinearFieldBuilder::Attribute("surface"); // TODO + +} // namespace katrin diff --git a/KEMField/Source/Bindings/Fields/Electric/src/KElectrostaticPotentialmapBuilder.cc b/KEMField/Source/Bindings/Fields/Electric/src/KElectrostaticPotentialmapBuilder.cc index 7ae4a98d3..35eceb2bf 100644 --- a/KEMField/Source/Bindings/Fields/Electric/src/KElectrostaticPotentialmapBuilder.cc +++ b/KEMField/Source/Bindings/Fields/Electric/src/KElectrostaticPotentialmapBuilder.cc @@ -35,6 +35,7 @@ STATICINT sKElectrostaticPotentialmapCalculatorStructure = KElectrostaticPotentialmapCalculatorBuilder::Attribute("directory") + KElectrostaticPotentialmapCalculatorBuilder::Attribute("file") + KElectrostaticPotentialmapCalculatorBuilder::Attribute("force_update") + + KElectrostaticPotentialmapCalculatorBuilder::Attribute("compute_field") + KElectrostaticPotentialmapCalculatorBuilder::Attribute("center") + KElectrostaticPotentialmapCalculatorBuilder::Attribute("length") + KElectrostaticPotentialmapCalculatorBuilder::Attribute("mirror_x") + diff --git a/KEMField/Source/Bindings/Fields/Magnetic/include/KMagnetostaticFieldmapBuilder.hh b/KEMField/Source/Bindings/Fields/Magnetic/include/KMagnetostaticFieldmapBuilder.hh index 662fcbb2a..eafa3551a 100644 --- a/KEMField/Source/Bindings/Fields/Magnetic/include/KMagnetostaticFieldmapBuilder.hh +++ b/KEMField/Source/Bindings/Fields/Magnetic/include/KMagnetostaticFieldmapBuilder.hh @@ -69,6 +69,10 @@ template<> inline bool KMagnetostaticFieldmapCalculatorBuilder::AddAttribute(KCo aContainer->CopyTo(fObject, &KEMField::KMagnetostaticFieldmapCalculator::SetForceUpdate); return true; } + if (aContainer->GetName() == "compute_gradient") { + aContainer->CopyTo(fObject, &KEMField::KMagnetostaticFieldmapCalculator::SetComputeGradient); + return true; + } if (aContainer->GetName() == "center") { KEMField::KEMStreamableThreeVector center; aContainer->CopyTo(center); diff --git a/KEMField/Source/Bindings/Fields/Magnetic/include/KRampedMagneticFieldBuilder.hh b/KEMField/Source/Bindings/Fields/Magnetic/include/KRampedMagneticFieldBuilder.hh index a1c38b532..1b1fadf7f 100644 --- a/KEMField/Source/Bindings/Fields/Magnetic/include/KRampedMagneticFieldBuilder.hh +++ b/KEMField/Source/Bindings/Fields/Magnetic/include/KRampedMagneticFieldBuilder.hh @@ -65,6 +65,9 @@ template<> inline bool KRampedMagneticFieldBuilder::AddAttribute(KContainer* aCo else if (aContainer->GetName() == "time_constant") { aContainer->CopyTo(fObject, &KEMField::KRampedMagneticField::SetTimeConstant); } + else if (aContainer->GetName() == "time_constant_2") { + aContainer->CopyTo(fObject, &KEMField::KRampedMagneticField::SetTimeConstant2); + } else if (aContainer->GetName() == "time_scaling") { aContainer->CopyTo(fObject, &KEMField::KRampedMagneticField::SetTimeScalingFactor); } diff --git a/KEMField/Source/Bindings/Fields/Magnetic/src/KMagnetostaticFieldmapBuilder.cc b/KEMField/Source/Bindings/Fields/Magnetic/src/KMagnetostaticFieldmapBuilder.cc index eb6489929..bc9e38dbd 100644 --- a/KEMField/Source/Bindings/Fields/Magnetic/src/KMagnetostaticFieldmapBuilder.cc +++ b/KEMField/Source/Bindings/Fields/Magnetic/src/KMagnetostaticFieldmapBuilder.cc @@ -34,6 +34,7 @@ STATICINT sKMagnetostaticFieldmapCalculatorStructure = KMagnetostaticFieldmapCalculatorBuilder::Attribute("directory") + KMagnetostaticFieldmapCalculatorBuilder::Attribute("file") + KMagnetostaticFieldmapCalculatorBuilder::Attribute("force_update") + + KMagnetostaticFieldmapCalculatorBuilder::Attribute("compute_gradient") + KMagnetostaticFieldmapCalculatorBuilder::Attribute("center") + KMagnetostaticFieldmapCalculatorBuilder::Attribute("length") + KMagnetostaticFieldmapCalculatorBuilder::Attribute("mirror_x") + diff --git a/KEMField/Source/FieldSolvers/ZonalHarmonic/Generator/include/KZonalHarmonicContainer.hh b/KEMField/Source/FieldSolvers/ZonalHarmonic/Generator/include/KZonalHarmonicContainer.hh index e23fd8948..f6e9c4659 100644 --- a/KEMField/Source/FieldSolvers/ZonalHarmonic/Generator/include/KZonalHarmonicContainer.hh +++ b/KEMField/Source/FieldSolvers/ZonalHarmonic/Generator/include/KZonalHarmonicContainer.hh @@ -43,6 +43,7 @@ template class KZonalHarmonicContainer { return fCoordinateSystem; } + SourcePointVector& GetCentralSourcePoints() { return fCentralSourcePoints; @@ -51,6 +52,7 @@ template class KZonalHarmonicContainer { return fCentralSourcePoints; } + SourcePointVector& GetRemoteSourcePoints() { return fRemoteSourcePoints; @@ -59,6 +61,11 @@ template class KZonalHarmonicContainer { return fRemoteSourcePoints; } + + std::set> CentralSourcePoints(); + std::set> RemoteSourcePoints(); + + ZonalHarmonicContainerVector& GetSubContainers() { return fSubContainers; @@ -335,6 +342,42 @@ template void KZonalHarmonicContainer::ComputeCoefficients(i } } +template std::set> KZonalHarmonicContainer::CentralSourcePoints() +{ + std::set> SPs; + + for (auto& sp : fCentralSourcePoints) { + auto z0 = sp->GetZ0() + fCoordinateSystem.GetOrigin().Z(); + auto rho = sp->GetRho(); + SPs.insert({z0, rho}); + } + + for (auto& subcontainer : fSubContainers) { + for (auto& sp : subcontainer->CentralSourcePoints()) { + SPs.insert(sp); + } + } + return SPs; +} + +template std::set> KZonalHarmonicContainer::RemoteSourcePoints() +{ + std::set> SPs; + + for (auto& sp : fRemoteSourcePoints) { + auto z0 = sp->GetZ0() + fCoordinateSystem.GetOrigin().Z(); + auto rho = sp->GetRho(); + SPs.insert({z0, rho}); + } + + for (auto& subcontainer : fSubContainers) { + for (auto& sp : subcontainer->RemoteSourcePoints()) { + SPs.insert(sp); + } + } + return SPs; +} + } // end namespace KEMField #endif /* KZONALHARMONICCONTAINER_DEF */ diff --git a/KEMField/Source/FieldSolvers/ZonalHarmonic/Generator/src/KZHCoefficientGeneratorElement.cc b/KEMField/Source/FieldSolvers/ZonalHarmonic/Generator/src/KZHCoefficientGeneratorElement.cc index 1f4609a53..98b4096c6 100644 --- a/KEMField/Source/FieldSolvers/ZonalHarmonic/Generator/src/KZHCoefficientGeneratorElement.cc +++ b/KEMField/Source/FieldSolvers/ZonalHarmonic/Generator/src/KZHCoefficientGeneratorElement.cc @@ -14,12 +14,13 @@ bool KZHCoefficientGeneratorElement::IsCoaxial(const KEMCoordinateSystem& coordi if (betweenOrigins.MagnitudeSquared() < coaxialityTolerance) { return true; } - else { - betweenOrigins = betweenOrigins.Unit(); - if (1. - fabs(GetCoordinateSystem().GetZAxis().Dot(betweenOrigins)) > coaxialityTolerance) - return false; - } + betweenOrigins = betweenOrigins.Unit(); + if (1. - fabs(GetCoordinateSystem().GetZAxis().Dot(betweenOrigins)) > coaxialityTolerance) + return false; + + if (coaxialityTolerance <= 1e-14) + return false; return true; } diff --git a/KEMField/Source/FieldSolvers/ZonalHarmonic/Solver/include/KElectromagnetZonalHarmonicFieldSolver.hh b/KEMField/Source/FieldSolvers/ZonalHarmonic/Solver/include/KElectromagnetZonalHarmonicFieldSolver.hh index d8354ff5e..d70933c8e 100644 --- a/KEMField/Source/FieldSolvers/ZonalHarmonic/Solver/include/KElectromagnetZonalHarmonicFieldSolver.hh +++ b/KEMField/Source/FieldSolvers/ZonalHarmonic/Solver/include/KElectromagnetZonalHarmonicFieldSolver.hh @@ -30,6 +30,9 @@ template<> class KZonalHarmonicFieldSolver : public KZonalH KGradient MagneticFieldGradient(const KPosition& P) const; std::pair MagneticFieldAndGradient(const KPosition& P) const; + bool CentralExpansion(const KPosition& P) const; + bool RemoteExpansion(const KPosition& P) const; + private: KZHLegendreCoefficients* fZHCoeffSingleton; diff --git a/KEMField/Source/FieldSolvers/ZonalHarmonic/Solver/include/KElectrostaticZonalHarmonicFieldSolver.hh b/KEMField/Source/FieldSolvers/ZonalHarmonic/Solver/include/KElectrostaticZonalHarmonicFieldSolver.hh index 69ae4e5d5..e881f407d 100644 --- a/KEMField/Source/FieldSolvers/ZonalHarmonic/Solver/include/KElectrostaticZonalHarmonicFieldSolver.hh +++ b/KEMField/Source/FieldSolvers/ZonalHarmonic/Solver/include/KElectrostaticZonalHarmonicFieldSolver.hh @@ -24,6 +24,9 @@ template<> class KZonalHarmonicFieldSolver : public KZonalH ~KZonalHarmonicFieldSolver() override {} + bool CentralExpansion(const KPosition& P) const; + bool RemoteExpansion(const KPosition& P) const; + double Potential(const KPosition& P) const; KThreeVector ElectricField(const KPosition& P) const; std::pair ElectricFieldAndPotential(const KPosition& P) const; diff --git a/KEMField/Source/FieldSolvers/ZonalHarmonic/Solver/src/KElectromagnetZonalHarmonicFieldSolver.cc b/KEMField/Source/FieldSolvers/ZonalHarmonic/Solver/src/KElectromagnetZonalHarmonicFieldSolver.cc index adce8b388..8cc1100f4 100644 --- a/KEMField/Source/FieldSolvers/ZonalHarmonic/Solver/src/KElectromagnetZonalHarmonicFieldSolver.cc +++ b/KEMField/Source/FieldSolvers/ZonalHarmonic/Solver/src/KElectromagnetZonalHarmonicFieldSolver.cc @@ -4,6 +4,37 @@ namespace KEMField { + +bool KZonalHarmonicFieldSolver::CentralExpansion(const KPosition& P) const +{ + KThreeVector localP = fContainer.GetCoordinateSystem().ToLocal(P); + + if (UseCentralExpansion(localP)) + return true; + + for (auto& sub : fSubsetFieldSolvers) { + if (sub->CentralExpansion(P)) + return true; + } + + return false; +} + +bool KZonalHarmonicFieldSolver::RemoteExpansion(const KPosition& P) const +{ + KThreeVector localP = fContainer.GetCoordinateSystem().ToLocal(P); + + if (UseRemoteExpansion(localP)) + return true; + + for (auto& sub : fSubsetFieldSolvers) { + if (sub->RemoteExpansion(P)) + return true; + } + + return false; +} + KThreeVector KZonalHarmonicFieldSolver::VectorPotential(const KPosition& P) const { KThreeVector localP = fContainer.GetCoordinateSystem().ToLocal(P); @@ -50,7 +81,6 @@ KThreeVector KZonalHarmonicFieldSolver::VectorPotential(con KThreeVector KZonalHarmonicFieldSolver::MagneticField(const KPosition& P) const { KThreeVector localP = fContainer.GetCoordinateSystem().ToLocal(P); - KThreeVector B; if (fCentralFirst) { @@ -969,8 +999,8 @@ bool KZonalHarmonicFieldSolver::CentralGradientExpansionNum for (unsigned int i = 0; i < 3; i++) { KThreeVector forward; KThreeVector backward; - retval *= CentralExpansionMagneticField(P + num_eps * coord_ax[i], forward); - retval *= CentralExpansionMagneticField(P - num_eps * coord_ax[i], backward); + retval = retval && CentralExpansionMagneticField(P + num_eps * coord_ax[i], forward); + retval = retval && CentralExpansionMagneticField(P - num_eps * coord_ax[i], backward); result[i] = (forward - backward) / (2.0 * num_eps); g(0, i) = result[i][0]; g(1, i) = result[i][1]; @@ -995,8 +1025,8 @@ bool KZonalHarmonicFieldSolver::RemoteGradientExpansionNume for (unsigned int i = 0; i < 3; i++) { KThreeVector forward; KThreeVector backward; - retval *= RemoteExpansionMagneticField(P + num_eps * coord_ax[i], forward); - retval *= RemoteExpansionMagneticField(P - num_eps * coord_ax[i], backward); + retval = retval && RemoteExpansionMagneticField(P + num_eps * coord_ax[i], forward); + retval = retval && RemoteExpansionMagneticField(P - num_eps * coord_ax[i], backward); result[i] = (forward - backward) / (2.0 * num_eps); g(0, i) = result[i][0]; g(1, i) = result[i][1]; diff --git a/KEMField/Source/FieldSolvers/ZonalHarmonic/Solver/src/KElectrostaticZonalHarmonicFieldSolver.cc b/KEMField/Source/FieldSolvers/ZonalHarmonic/Solver/src/KElectrostaticZonalHarmonicFieldSolver.cc index af9b8d7ce..4a58b4a3c 100644 --- a/KEMField/Source/FieldSolvers/ZonalHarmonic/Solver/src/KElectrostaticZonalHarmonicFieldSolver.cc +++ b/KEMField/Source/FieldSolvers/ZonalHarmonic/Solver/src/KElectrostaticZonalHarmonicFieldSolver.cc @@ -6,6 +6,32 @@ namespace KEMField { +bool KZonalHarmonicFieldSolver::CentralExpansion(const KPosition& P) const +{ + if (UseCentralExpansion(P)) + return true; + + for (auto& sub : fSubsetFieldSolvers) { + if (sub->CentralExpansion(P)) + return true; + } + + return false; +} + +bool KZonalHarmonicFieldSolver::RemoteExpansion(const KPosition& P) const +{ + if (UseRemoteExpansion(P)) + return true; + + for (auto& sub : fSubsetFieldSolvers) { + if (sub->UseRemoteExpansion(P)) + return true; + } + + return false; +} + double KZonalHarmonicFieldSolver::Potential(const KPosition& P) const { double phi = 0; diff --git a/KEMField/Source/Interface/FieldSolvers/Electric/include/KElectricZHFieldSolver.hh b/KEMField/Source/Interface/FieldSolvers/Electric/include/KElectricZHFieldSolver.hh index 8353d0551..f51931aaf 100644 --- a/KEMField/Source/Interface/FieldSolvers/Electric/include/KElectricZHFieldSolver.hh +++ b/KEMField/Source/Interface/FieldSolvers/Electric/include/KElectricZHFieldSolver.hh @@ -20,6 +20,8 @@ namespace KEMField class KElectricZHFieldSolver : public KElectricFieldSolver { public: + typedef std::vector SourcePointVector; + KElectricZHFieldSolver(); ~KElectricZHFieldSolver() override; @@ -36,6 +38,9 @@ class KElectricZHFieldSolver : public KElectricFieldSolver return fParameters; } + std::set> CentralSourcePoints(); + std::set> RemoteSourcePoints(); + private: void InitializeCore(KSurfaceContainer& container) override; diff --git a/KEMField/Source/Interface/FieldSolvers/Electric/src/KElectricZHFieldSolver.cc b/KEMField/Source/Interface/FieldSolvers/Electric/src/KElectricZHFieldSolver.cc index 1bb767476..a947ef877 100644 --- a/KEMField/Source/Interface/FieldSolvers/Electric/src/KElectricZHFieldSolver.cc +++ b/KEMField/Source/Interface/FieldSolvers/Electric/src/KElectricZHFieldSolver.cc @@ -121,12 +121,22 @@ std::pair KElectricZHFieldSolver::ElectricFieldAndPotentia bool KElectricZHFieldSolver::UseCentralExpansion(const KPosition& P) { - return fZonalHarmonicFieldSolver->UseCentralExpansion(P); + return fZonalHarmonicFieldSolver->CentralExpansion(P); } bool KElectricZHFieldSolver::UseRemoteExpansion(const KPosition& P) { - return fZonalHarmonicFieldSolver->UseRemoteExpansion(P); + return fZonalHarmonicFieldSolver->RemoteExpansion(P); +} + +std::set> KElectricZHFieldSolver::CentralSourcePoints() +{ + return fZHContainer->CentralSourcePoints(); +} + +std::set> KElectricZHFieldSolver::RemoteSourcePoints() +{ + return fZHContainer->RemoteSourcePoints(); } diff --git a/KEMField/Source/Interface/FieldSolvers/Magnetic/include/KZonalHarmonicMagnetostaticFieldSolver.hh b/KEMField/Source/Interface/FieldSolvers/Magnetic/include/KZonalHarmonicMagnetostaticFieldSolver.hh index 0b4ac5a53..20f0ea1c4 100644 --- a/KEMField/Source/Interface/FieldSolvers/Magnetic/include/KZonalHarmonicMagnetostaticFieldSolver.hh +++ b/KEMField/Source/Interface/FieldSolvers/Magnetic/include/KZonalHarmonicMagnetostaticFieldSolver.hh @@ -19,9 +19,17 @@ namespace KEMField class KZonalHarmonicMagnetostaticFieldSolver : public KMagneticFieldSolver { public: + typedef std::vector SourcePointVector; + KZonalHarmonicMagnetostaticFieldSolver(); ~KZonalHarmonicMagnetostaticFieldSolver() override; + bool UseCentralExpansion(const KPosition& P); + bool UseRemoteExpansion(const KPosition& P); + + std::set> CentralSourcePoints(); + std::set> RemoteSourcePoints(); + void InitializeCore(KElectromagnetContainer& container) override; KThreeVector MagneticPotentialCore(const KPosition& P) const override; @@ -34,6 +42,7 @@ class KZonalHarmonicMagnetostaticFieldSolver : public KMagneticFieldSolver return fParameters; } + private: KElectromagnetIntegrator fIntegrator; KZonalHarmonicContainer* fZHContainer; diff --git a/KEMField/Source/Interface/FieldSolvers/Magnetic/src/KZonalHarmonicMagnetostaticFieldSolver.cc b/KEMField/Source/Interface/FieldSolvers/Magnetic/src/KZonalHarmonicMagnetostaticFieldSolver.cc index 553a6dcc1..e87322148 100644 --- a/KEMField/Source/Interface/FieldSolvers/Magnetic/src/KZonalHarmonicMagnetostaticFieldSolver.cc +++ b/KEMField/Source/Interface/FieldSolvers/Magnetic/src/KZonalHarmonicMagnetostaticFieldSolver.cc @@ -97,4 +97,25 @@ KZonalHarmonicMagnetostaticFieldSolver::MagneticFieldAndGradientCore(const KPosi return fZonalHarmonicFieldSolver->MagneticFieldAndGradient(P); } +bool KZonalHarmonicMagnetostaticFieldSolver::UseCentralExpansion(const KPosition& P) +{ + return fZonalHarmonicFieldSolver->CentralExpansion(P); +} + +bool KZonalHarmonicMagnetostaticFieldSolver::UseRemoteExpansion(const KPosition& P) +{ + return fZonalHarmonicFieldSolver->RemoteExpansion(P); +} + +std::set> KZonalHarmonicMagnetostaticFieldSolver::CentralSourcePoints() +{ + return fZHContainer->CentralSourcePoints(); +} + +std::set> KZonalHarmonicMagnetostaticFieldSolver::RemoteSourcePoints() +{ + return fZHContainer->RemoteSourcePoints(); +} + + } /* namespace KEMField */ diff --git a/KEMField/Source/Interface/Fields/Electric/CMakeLists.txt b/KEMField/Source/Interface/Fields/Electric/CMakeLists.txt index ce2b0d1b4..3a2c95cc1 100644 --- a/KEMField/Source/Interface/Fields/Electric/CMakeLists.txt +++ b/KEMField/Source/Interface/Fields/Electric/CMakeLists.txt @@ -7,6 +7,7 @@ set( FIELDS_ELECTRIC_HEADER_BASENAMES KElectricField.hh KElectricQuadrupoleField.hh KElectrostaticConstantField.hh + KElectrostaticLinearField.hh KElectrostaticField.hh KElectrostaticBoundaryField.hh KElectrostaticBoundaryFieldTimer.hh @@ -27,6 +28,7 @@ set( FIELDS_ELECTRIC_SOURCE_BASENAMES KBoundaryElementInfoDisplay.cc KElectricQuadrupoleField.cc KElectrostaticConstantField.cc + KElectrostaticLinearField.cc KElectrostaticBoundaryField.cc KElectrostaticBoundaryFieldTimer.cc KInducedAzimuthalElectricField.cc diff --git a/KEMField/Source/Interface/Fields/Electric/include/KElectrostaticLinearField.hh b/KEMField/Source/Interface/Fields/Electric/include/KElectrostaticLinearField.hh new file mode 100644 index 000000000..dd740c2e0 --- /dev/null +++ b/KEMField/Source/Interface/Fields/Electric/include/KElectrostaticLinearField.hh @@ -0,0 +1,50 @@ +#ifndef KELECTROSTATICLINEARFIELD_DEF +#define KELECTROSTATICLINEARFIELD_DEF + +#include "KElectrostaticField.hh" +#include "KGCore.hh" + +namespace KEMField +{ + +class KElectrostaticLinearField : public KElectrostaticField +{ + public: + KElectrostaticLinearField(); + ~KElectrostaticLinearField() override{}; + + static std::string Name() + { + return "ElectrostaticLinearFieldSolver"; + } + + private: + double PotentialCore(const KPosition& aSamplePoint) const override; + KThreeVector ElectricFieldCore(const KPosition& aSamplePoint) const override; + + public: + void SetPotential1(double aPotential); + double GetPotential1() const; + + void SetPotential2(double aPotential); + double GetPotential2() const; + + void SetZ1(double aPosition); + double GetZ1() const; + + void SetZ2(double aPosition); + double GetZ2() const; + + void SetSurface(KGeoBag::KGSurface* aSurface); + const KGeoBag::KGSurface* GetSurface() const; + + + protected: + double fU1, fU2; + double fZ1, fZ2; + KGeoBag::KGSurface* fSurface; +}; + +} // namespace KEMField + +#endif /* KELECTROSTATICLINEARFIELD_DEF */ diff --git a/KEMField/Source/Interface/Fields/Electric/src/KElectrostaticLinearField.cc b/KEMField/Source/Interface/Fields/Electric/src/KElectrostaticLinearField.cc new file mode 100644 index 000000000..1ca0c1073 --- /dev/null +++ b/KEMField/Source/Interface/Fields/Electric/src/KElectrostaticLinearField.cc @@ -0,0 +1,89 @@ +#include "KElectrostaticLinearField.hh" + +#include "KGVolume.hh" + +#include + +namespace KEMField +{ + +KElectrostaticLinearField::KElectrostaticLinearField() : + KElectrostaticField(), + fU1(0), + fU2(0), + fZ1(0), + fZ2(0), + fSurface(nullptr) +{} + +double KElectrostaticLinearField::PotentialCore(const KPosition& aSamplePoint) const +{ + // TODO: coordinate transform + + if ((aSamplePoint.Z() < fZ1) || (aSamplePoint.Z()) > fZ2) + return 0; + + auto E = (fU2 - fU1) / (fZ2 - fZ1); + + return (aSamplePoint.Z() - fZ1) * E + fU1; +} + +KThreeVector KElectrostaticLinearField::ElectricFieldCore(const KPosition& aSamplePoint) const +{ + // TODO: coordinate transform + + if ((aSamplePoint.Z() < fZ1) || (aSamplePoint.Z()) > fZ2) + return KThreeVector::sZero; + + auto E = (fU2 - fU1) / (fZ2 - fZ1); + + return KThreeVector(0, 0, E); +} + +void KElectrostaticLinearField::SetPotential1(double aPotential) +{ + fU1 = aPotential; +} +double KElectrostaticLinearField::GetPotential1() const +{ + return fU1; +} + +void KElectrostaticLinearField::SetPotential2(double aPotential) +{ + fU2 = aPotential; +} +double KElectrostaticLinearField::GetPotential2() const +{ + return fU2; +} + +void KElectrostaticLinearField::SetZ1(double aPosition) +{ + fZ1 = aPosition; +} +double KElectrostaticLinearField::GetZ1() const +{ + return fZ1; +} + +void KElectrostaticLinearField::SetZ2(double aPosition) +{ + fZ2 = aPosition; +} +double KElectrostaticLinearField::GetZ2() const +{ + return fZ2; +} + + +void KElectrostaticLinearField::SetSurface(KGeoBag::KGSurface* aSurface) +{ + fSurface = aSurface; +} +const KGeoBag::KGSurface* KElectrostaticLinearField::GetSurface() const +{ + return fSurface; +} + +} // namespace KEMField diff --git a/KEMField/Source/Interface/KGeoBag/include/KGBEMConverter.hh b/KEMField/Source/Interface/KGeoBag/include/KGBEMConverter.hh index 92b080093..f56610b52 100644 --- a/KEMField/Source/Interface/KGeoBag/include/KGBEMConverter.hh +++ b/KEMField/Source/Interface/KGeoBag/include/KGBEMConverter.hh @@ -312,7 +312,7 @@ class KGBEMMeshConverter : public KGDualHierarchy @@ -327,7 +327,7 @@ class KGBEMAxialMeshConverter : public KGDualHierarchy @@ -342,7 +342,7 @@ class KGBEMDiscreteRotationalMeshConverter : public KGDualHierarchyAsExtension()); + if (!Add(aSurface->AsExtension())) + coremsg(eWarning) << "not adding surface <" << aSurface->GetPath() << "> since it is not coaxial" << eom; return; } void KGBEMAxialMeshConverter::DispatchSpace(KGSpace* aSpace) { - Add(aSpace->AsExtension()); + if (!Add(aSpace->AsExtension())) + coremsg(eWarning) << "not adding space <" << aSpace->GetPath() << "> since it is not coaxial" << eom; return; } -void KGBEMAxialMeshConverter::Add(KGAxialMeshData* aData) +bool KGBEMAxialMeshConverter::Add(KGAxialMeshData* aData) { KGAxialMeshElement* tAxialMeshElement; KGAxialMeshLoop* tAxialMeshLoop; @@ -309,7 +311,7 @@ void KGBEMAxialMeshConverter::Add(KGAxialMeshData* aData) //cout << "...current origin is <" << fCurrentOrigin << ">" << endl; //cout << "...current z axis is <" << fCurrentZAxis << ">" << endl; //cout << "...axes do not match!" << endl; - return; + return false; } for (auto tElementIt = aData->Elements()->begin(); tElementIt != aData->Elements()->end(); tElementIt++) { @@ -336,7 +338,7 @@ void KGBEMAxialMeshConverter::Add(KGAxialMeshData* aData) //cout << "...added <" << fConicSections.size() + fRings.size() << "> components." << endl; } - return; + return true; } KGBEMDiscreteRotationalMeshConverter::KGBEMDiscreteRotationalMeshConverter() {} @@ -348,16 +350,18 @@ KGBEMDiscreteRotationalMeshConverter::~KGBEMDiscreteRotationalMeshConverter() {} void KGBEMDiscreteRotationalMeshConverter::DispatchSurface(KGSurface* aSurface) { - Add(aSurface->AsExtension()); + if (!Add(aSurface->AsExtension())) + coremsg(eWarning) << "not adding surface <" << aSurface->GetPath() << "> since it is not coaxial" << eom; return; } void KGBEMDiscreteRotationalMeshConverter::DispatchSpace(KGSpace* aSpace) { - Add(aSpace->AsExtension()); + if (!Add(aSpace->AsExtension())) + coremsg(eWarning) << "not adding space <" << aSpace->GetPath() << "> since it is not coaxial" << eom; return; } -void KGBEMDiscreteRotationalMeshConverter::Add(KGDiscreteRotationalMeshData* aData) +bool KGBEMDiscreteRotationalMeshConverter::Add(KGDiscreteRotationalMeshData* aData) { KGDiscreteRotationalMeshElement* tMeshElement; KGDiscreteRotationalMeshRectangle* tMeshRectangle; @@ -377,7 +381,7 @@ void KGBEMDiscreteRotationalMeshConverter::Add(KGDiscreteRotationalMeshData* aDa if (fAxis.EqualTo(fCurrentAxis) == false) { // improve the hell out of this //cout << "...axes do not match!" << endl; - return; + return false; } tCenter.SetComponents(fAxis.GetCenter().X(), fAxis.GetCenter().Y(), fAxis.GetCenter().Z()); @@ -425,6 +429,6 @@ void KGBEMDiscreteRotationalMeshConverter::Add(KGDiscreteRotationalMeshData* aDa } } - return; + return true; } } // namespace KGeoBag diff --git a/KEMField/Source/Interface/KGeoBag/src/KGElectromagnetConverter.cc b/KEMField/Source/Interface/KGeoBag/src/KGElectromagnetConverter.cc index 5849a9591..f7774675a 100644 --- a/KEMField/Source/Interface/KGeoBag/src/KGElectromagnetConverter.cc +++ b/KEMField/Source/Interface/KGeoBag/src/KGElectromagnetConverter.cc @@ -5,6 +5,9 @@ #include "KGCylinderSurface.hh" #include "KGRodSpace.hh" +/// FIXME: If defined, print magfield3 lines of converted magnet geometry. This should be a runtime option. +//#define PRINT_MAGFIELD3 + namespace KGeoBag { @@ -92,6 +95,10 @@ void KGElectromagnetConverter::VisitSpace(KGSpace* aSpace) { Clear(); +#ifdef PRINT_MAGFIELD3 + std::cout << "# " << aSpace->GetPath() << std::endl; +#endif + fCurrentOrigin = aSpace->GetOrigin(); fCurrentXAxis = aSpace->GetXAxis(); fCurrentYAxis = aSpace->GetYAxis(); @@ -114,6 +121,10 @@ void KGElectromagnetConverter::VisitSurface(KGSurface* aSurface) void KGElectromagnetConverter::VisitExtendedSpace(KGExtendedSpace* electromagnetSpace) { +#ifdef PRINT_MAGFIELD3 + //std::cout << "# " << electromagnetSpace->GetName() << std::endl; +#endif + fCurrentElectromagnetSpace = electromagnetSpace; } @@ -181,6 +192,16 @@ void KGElectromagnetConverter::VisitCylinderTubeSpace(KGCylinderTubeSpace* cylin GlobalToInternalVector(fCurrentYAxis), GlobalToInternalVector(fCurrentZAxis)); fElectromagnetContainer->push_back(coil); + +#ifdef PRINT_MAGFIELD3 + // do not use coil->GetP0|P1() because it is defined as (r,0,z) + auto p0 = coil->GetCoordinateSystem().ToGlobal(KPosition(0, 0, tZMin)); + auto p1 = coil->GetCoordinateSystem().ToGlobal(KPosition(0, 0, tZMax)); + + std::cout << " " << coil->GetCurrentDensity() << " " << p0.X() << " " << p0.Y() << " " << p0.Z() << " " + << p1.X() << " " << p1.Y() << " " << p1.Z() << " " << coil->GetR0() << " " << coil->GetR1() << " " + << coil->GetIntegrationScale() << std::endl; +#endif } } diff --git a/KEMField/Source/Interface/KGeoBag/src/KGElectrostaticBoundaryField.cc b/KEMField/Source/Interface/KGeoBag/src/KGElectrostaticBoundaryField.cc index 1ef57b897..ec654590e 100644 --- a/KEMField/Source/Interface/KGeoBag/src/KGElectrostaticBoundaryField.cc +++ b/KEMField/Source/Interface/KGeoBag/src/KGElectrostaticBoundaryField.cc @@ -128,10 +128,9 @@ void KGElectrostaticBoundaryField::ConfigureSurfaceContainer() } if (container->empty()) { - // TODO find alternative to KNamed::GetName which is not available here - cout << "ERROR:" - << "electrostatic field solver <" /*<< GetName()*/ << "> has zero surface elements" << endl; - std::exit(-1); + cout << "WARNING:" + << "electrostatic field solver <" << GetName() << "> has zero surface elements" << endl; + //std::exit(-1); } SetContainer(container); diff --git a/KEMField/Source/Interface/KGeoBag/src/KGStaticElectromagnetField.cc b/KEMField/Source/Interface/KGeoBag/src/KGStaticElectromagnetField.cc index a640087ca..decf227e3 100644 --- a/KEMField/Source/Interface/KGeoBag/src/KGStaticElectromagnetField.cc +++ b/KEMField/Source/Interface/KGeoBag/src/KGStaticElectromagnetField.cc @@ -93,9 +93,9 @@ void KGStaticElectromagnetField::ConfigureSurfaceContainer() } if (container->empty()) { - cout << "ERROR:" - << "electromagnet field solver <" /*<< GetName()*/ << "> has zero surface or space elements" << endl; - std::exit(-1); + cout << "WARNING:" + << "electromagnet field solver <" << GetName() << "> has zero surface or space elements" << endl; + //std::exit(-1); } SetContainer(container); diff --git a/KEMField/Source/Plugins/VTKPart2/include/KElectrostaticPotentialmap.hh b/KEMField/Source/Plugins/VTKPart2/include/KElectrostaticPotentialmap.hh index 5ec2eb967..5110a0f50 100644 --- a/KEMField/Source/Plugins/VTKPart2/include/KElectrostaticPotentialmap.hh +++ b/KEMField/Source/Plugins/VTKPart2/include/KElectrostaticPotentialmap.hh @@ -225,6 +225,10 @@ class KElectrostaticPotentialmapCalculator { fForceUpdate = aFlag; } + void SetComputeField(bool aFlag) + { + fComputeField = aFlag; + } void SetCenter(KPosition aCenter) { fCenter = aCenter; @@ -296,6 +300,7 @@ class KElectrostaticPotentialmapCalculator std::string fFile; std::string fName; bool fForceUpdate; + bool fComputeField; KPosition fCenter; KThreeVector fLength; bool fMirrorX, fMirrorY, fMirrorZ; diff --git a/KEMField/Source/Plugins/VTKPart2/include/KMagnetostaticFieldmap.hh b/KEMField/Source/Plugins/VTKPart2/include/KMagnetostaticFieldmap.hh index 641dfa512..ef5f0f072 100644 --- a/KEMField/Source/Plugins/VTKPart2/include/KMagnetostaticFieldmap.hh +++ b/KEMField/Source/Plugins/VTKPart2/include/KMagnetostaticFieldmap.hh @@ -228,6 +228,10 @@ class KMagnetostaticFieldmapCalculator { fForceUpdate = aFlag; } + void SetComputeGradient(bool aFlag) + { + fComputeGradient = aFlag; + } void SetCenter(KPosition aCenter) { fCenter = aCenter; @@ -299,6 +303,7 @@ class KMagnetostaticFieldmapCalculator std::string fFile; std::string fName; bool fForceUpdate; + bool fComputeGradient; KPosition fCenter; KThreeVector fLength; bool fMirrorX, fMirrorY, fMirrorZ; @@ -309,6 +314,7 @@ class KMagnetostaticFieldmapCalculator vtkSmartPointer fGrid; vtkSmartPointer fValidityData; vtkSmartPointer fFieldData; + vtkSmartPointer fGradientData; }; } /* namespace KEMField */ diff --git a/KEMField/Source/Plugins/VTKPart2/src/KElectrostaticPotentialmap.cc b/KEMField/Source/Plugins/VTKPart2/src/KElectrostaticPotentialmap.cc index ebdb6d7e5..fd35ee9d8 100644 --- a/KEMField/Source/Plugins/VTKPart2/src/KElectrostaticPotentialmap.cc +++ b/KEMField/Source/Plugins/VTKPart2/src/KElectrostaticPotentialmap.cc @@ -11,6 +11,7 @@ #include "KEMFileInterface.hh" #include "KEMSimpleException.hh" #include "KFile.h" +#include "KGslErrorHandler.h" #include #include @@ -388,6 +389,8 @@ KElectrostaticPotentialmapCalculator::KElectrostaticPotentialmapCalculator() : fOutputFilename(""), fDirectory(SCRATCH_DEFAULT_DIR), fFile(""), + fForceUpdate(false), + fComputeField(false), fCenter(), fLength(), fMirrorX(false), @@ -510,11 +513,13 @@ void KElectrostaticPotentialmapCalculator::Prepare() fPotentialData->SetNumberOfTuples(tNumPoints); fGrid->GetPointData()->AddArray(fPotentialData); - fFieldData = vtkSmartPointer::New(); - fFieldData->SetName("electric field"); - fFieldData->SetNumberOfComponents(3); // vector data - fFieldData->SetNumberOfTuples(tNumPoints); - fGrid->GetPointData()->AddArray(fFieldData); + if (fComputeField) { + fFieldData = vtkSmartPointer::New(); + fFieldData->SetName("electric field"); + fFieldData->SetNumberOfComponents(3); // vector data + fFieldData->SetNumberOfTuples(tNumPoints); + fGrid->GetPointData()->AddArray(fFieldData); + } } void KElectrostaticPotentialmapCalculator::Execute() @@ -577,8 +582,13 @@ void KElectrostaticPotentialmapCalculator::Execute() if (!tHasValue) { tPotential = 0.; - for (auto& it : fElectricFields) - tPotential += it.second->Potential(KThreeVector(tPoint)); + try { + for (auto& it : fElectricFields) + tPotential += it.second->Potential(KThreeVector(tPoint)); + } + catch (katrin::KGslException& e) { + continue; + } } fPotentialData->SetTuple1(i, tPotential); @@ -592,65 +602,76 @@ void KElectrostaticPotentialmapCalculator::Execute() << ", time per potential evaluation = " << tTimeSpent / (double) (tNumPoints) << ")" << endl; - cout << "computing electric field at " << tNumPoints << " grid points" << endl; - - //evaluate field - tClockStart = clock(); - for (unsigned int i = 0; i < tNumPoints; i++) { - if (i % 10 == 0) { - int progress = 50 * (float) i / (float) (tNumPoints - 1); - std::cout << "\r "; - for (int j = 0; j < 50; j++) - std::cout << (j <= progress ? "#" : "."); - std::cout << " [" << 2 * progress << "%]" << std::flush; - } - - double tPoint[3]; - fGrid->GetPoint(i, tPoint); - - if (!CheckPosition(KPosition(tPoint))) - continue; + if (fComputeField) { + cout << "computing electric field at " << tNumPoints << " grid points" << endl; - bool tHasValue = false; - KThreeVector tField; + //evaluate field + tClockStart = clock(); + for (unsigned int i = 0; i < tNumPoints; i++) { + if (i % 10 == 0) { + int progress = 50 * (float) i / (float) (tNumPoints - 1); + std::cout << "\r "; + for (int j = 0; j < 50; j++) + std::cout << (j <= progress ? "#" : "."); + std::cout << " [" << 2 * progress << "%]" << std::flush; + } - if (fMirrorX || fMirrorY || fMirrorZ) { - double tMirrorPoint[3]; - tMirrorPoint[0] = tPoint[0]; - tMirrorPoint[1] = tPoint[1]; - tMirrorPoint[2] = tPoint[2]; - if (fMirrorX && (tPoint[0] > fCenter.X())) - tMirrorPoint[0] = 2. * fCenter.X() - tPoint[0]; - if (fMirrorY && (tPoint[1] > fCenter.Y())) - tMirrorPoint[1] = 2. * fCenter.Y() - tPoint[1]; - if (fMirrorZ && (tPoint[2] > fCenter.Z())) - tMirrorPoint[2] = 2. * fCenter.Z() - tPoint[2]; + double tPoint[3]; + fGrid->GetPoint(i, tPoint); + + if (!CheckPosition(KPosition(tPoint))) + continue; + + bool tHasValue = false; + KThreeVector tField; + + if (fMirrorX || fMirrorY || fMirrorZ) { + double tMirrorPoint[3]; + tMirrorPoint[0] = tPoint[0]; + tMirrorPoint[1] = tPoint[1]; + tMirrorPoint[2] = tPoint[2]; + if (fMirrorX && (tPoint[0] > fCenter.X())) + tMirrorPoint[0] = 2. * fCenter.X() - tPoint[0]; + if (fMirrorY && (tPoint[1] > fCenter.Y())) + tMirrorPoint[1] = 2. * fCenter.Y() - tPoint[1]; + if (fMirrorZ && (tPoint[2] > fCenter.Z())) + tMirrorPoint[2] = 2. * fCenter.Z() - tPoint[2]; + + if ((tMirrorPoint[0] != tPoint[0]) || (tMirrorPoint[1] != tPoint[1]) || + (tMirrorPoint[2] != tPoint[2])) { + unsigned int j = fGrid->FindPoint(tMirrorPoint); + if (fValidityData->GetTuple1(j) >= 2) // 2 = field valid + { + tField.SetComponents(fFieldData->GetTuple3(j)); + tHasValue = true; + } + } + } - if ((tMirrorPoint[0] != tPoint[0]) || (tMirrorPoint[1] != tPoint[1]) || (tMirrorPoint[2] != tPoint[2])) { - unsigned int j = fGrid->FindPoint(tMirrorPoint); - if (fValidityData->GetTuple1(j) >= 2) // 2 = field valid - { - tField.SetComponents(fFieldData->GetTuple3(j)); - tHasValue = true; + if (!tHasValue) { + tField = KThreeVector::sZero; + try { + for (auto& it : fElectricFields) + tField += it.second->ElectricField(KPosition(tPoint)); + } + catch (katrin::KGslException& e) { + continue; } } - } - if (!tHasValue) { - tField = KThreeVector::sZero; - for (auto& it : fElectricFields) - tField += it.second->ElectricField(KPosition(tPoint)); + fFieldData->SetTuple3(i, tField[0], tField[1], tField[2]); + fValidityData->SetTuple1(i, 2); } + std::cout << std::endl; + tClockEnd = clock(); - fFieldData->SetTuple3(i, tField[0], tField[1], tField[2]); - fValidityData->SetTuple1(i, 2); + tTimeSpent = ((double) (tClockEnd - tClockStart)) / CLOCKS_PER_SEC; // time in seconds + cout << "finished computing field map (total time spent = " << tTimeSpent + << ", time per field evaluation = " << tTimeSpent / (double) (tNumPoints) << ")" << endl; + } + else { + cout << "not computing electric field" << endl; } - std::cout << std::endl; - tClockEnd = clock(); - - tTimeSpent = ((double) (tClockEnd - tClockStart)) / CLOCKS_PER_SEC; // time in seconds - cout << "finished computing potential map (total time spent = " << tTimeSpent - << ", time per field evaluation = " << tTimeSpent / (double) (tNumPoints) << ")" << endl; } void KElectrostaticPotentialmapCalculator::Finish() diff --git a/KEMField/Source/Plugins/VTKPart2/src/KMagnetostaticFieldmap.cc b/KEMField/Source/Plugins/VTKPart2/src/KMagnetostaticFieldmap.cc index 7cbcf12ac..7c297d549 100644 --- a/KEMField/Source/Plugins/VTKPart2/src/KMagnetostaticFieldmap.cc +++ b/KEMField/Source/Plugins/VTKPart2/src/KMagnetostaticFieldmap.cc @@ -11,6 +11,7 @@ #include "KEMFileInterface.hh" #include "KEMSimpleException.hh" #include "KFile.h" +#include "KGslErrorHandler.h" #include #include @@ -394,6 +395,8 @@ KMagnetostaticFieldmapCalculator::KMagnetostaticFieldmapCalculator() : fOutputFilename(""), fDirectory(SCRATCH_DEFAULT_DIR), fFile(""), + fForceUpdate(false), + fComputeGradient(false), fCenter(), fLength(), fMirrorX(false), @@ -514,6 +517,18 @@ void KMagnetostaticFieldmapCalculator::Prepare() fFieldData->SetNumberOfComponents(3); // vector data fFieldData->SetNumberOfTuples(tNumPoints); fGrid->GetPointData()->AddArray(fFieldData); + + if (fComputeGradient) { + if (fMagneticFields.size() > 1) { + cout << "computing magnetic gradient of multiple fields, results are probably incorrect" << endl; + } + + fGradientData = vtkSmartPointer::New(); + fGradientData->SetName("magnetic gradient"); + fGradientData->SetNumberOfComponents(9); // tensor data + fGradientData->SetNumberOfTuples(tNumPoints); + fGrid->GetPointData()->AddArray(fGradientData); + } } void KMagnetostaticFieldmapCalculator::Execute() @@ -576,8 +591,13 @@ void KMagnetostaticFieldmapCalculator::Execute() if (!tHasValue) { tField = KThreeVector::sZero; - for (auto& it : fMagneticFields) - tField += it.second->MagneticField(KPosition(tPoint)); + try { + for (auto& it : fMagneticFields) + tField += it.second->MagneticField(KPosition(tPoint)); + } + catch (katrin::KGslException& e) { + tField = KThreeVector::sInvalid; + } } fFieldData->SetTuple3(i, tField[0], tField[1], tField[2]); @@ -589,6 +609,89 @@ void KMagnetostaticFieldmapCalculator::Execute() tTimeSpent = ((double) (tClockEnd - tClockStart)) / CLOCKS_PER_SEC; // time in seconds cout << "finished computing field map (total time spent = " << tTimeSpent << ", time per field evaluation = " << tTimeSpent / (double) (tNumPoints) << ")" << endl; + + + if (fComputeGradient) { + cout << "computing magnetic gradient at " << tNumPoints << " grid points" << endl; + + //evaluate field + tClockStart = clock(); + for (unsigned int i = 0; i < tNumPoints; i++) { + if (i % 10 == 0) { + int progress = 50 * (float) i / (float) (tNumPoints - 1); + std::cout << "\r "; + for (int j = 0; j < 50; j++) + std::cout << (j <= progress ? "#" : "."); + std::cout << " [" << 2 * progress << "%]" << std::flush; + } + + double tPoint[3]; + fGrid->GetPoint(i, tPoint); + + if (!CheckPosition(KPosition(tPoint))) + continue; + + bool tHasValue = false; + KThreeMatrix tGradient; + + if (fMirrorX || fMirrorY || fMirrorZ) { + double tMirrorPoint[3]; + tMirrorPoint[0] = tPoint[0]; + tMirrorPoint[1] = tPoint[1]; + tMirrorPoint[2] = tPoint[2]; + if (fMirrorX && (tPoint[0] > fCenter.X())) + tMirrorPoint[0] = 2. * fCenter.X() - tPoint[0]; + if (fMirrorY && (tPoint[1] > fCenter.Y())) + tMirrorPoint[1] = 2. * fCenter.Y() - tPoint[1]; + if (fMirrorZ && (tPoint[2] > fCenter.Z())) + tMirrorPoint[2] = 2. * fCenter.Z() - tPoint[2]; + + if ((tMirrorPoint[0] != tPoint[0]) || (tMirrorPoint[1] != tPoint[1]) || + (tMirrorPoint[2] != tPoint[2])) { + unsigned int j = fGrid->FindPoint(tMirrorPoint); + if (fValidityData->GetTuple1(j) >= 3) // 3 = gradient valid + { + tGradient.SetComponents(fGradientData->GetTuple9(j)); + tHasValue = true; + } + } + } + + if (!tHasValue) { + tGradient = KThreeMatrix::sZero; + try { + /// FIXME: summing up gradients from multiple fields doesn't make much sense, + /// Needs to be handled in a better way, e.g. calculate gradient from field map directly. + for (auto& it : fMagneticFields) + tGradient += it.second->MagneticGradient(KPosition(tPoint)); + } + catch (katrin::KGslException& e) { + tGradient = KThreeMatrix::sInvalid; + } + } + + fGradientData->SetTuple9(i, + tGradient[0], + tGradient[1], + tGradient[2], + tGradient[3], + tGradient[4], + tGradient[5], + tGradient[6], + tGradient[7], + tGradient[8]); + fValidityData->SetTuple1(i, 3); + } + std::cout << std::endl; + tClockEnd = clock(); + + tTimeSpent = ((double) (tClockEnd - tClockStart)) / CLOCKS_PER_SEC; // time in seconds + cout << "finished computing gradient map (total time spent = " << tTimeSpent + << ", time per gradient evaluation = " << tTimeSpent / (double) (tNumPoints) << ")" << endl; + } + else { + cout << "not computing magnetic gradient" << endl; + } } void KMagnetostaticFieldmapCalculator::Finish() diff --git a/KGeoBag/CMakeLists.txt b/KGeoBag/CMakeLists.txt index 66e3f8b20..4dbfdbd77 100644 --- a/KGeoBag/CMakeLists.txt +++ b/KGeoBag/CMakeLists.txt @@ -3,7 +3,7 @@ cmake_minimum_required( VERSION ${CMAKE_MINIMUM_VERSION} ) # KGeoBag version set( MODULE_VERSION_MAJOR 3 ) set( MODULE_VERSION_MINOR 7 ) -set( MODULE_VERSION_PATCH 0 ) +set( MODULE_VERSION_PATCH 3 ) set( MODULE_VERSION "${MODULE_VERSION_MAJOR}.${MODULE_VERSION_MINOR}.${MODULE_VERSION_PATCH}" ) #project( KGeoBag VERSION ${MODULE_VERSION}) diff --git a/KGeoBag/Source/Application/Source/KGeoBag.cc b/KGeoBag/Source/Application/Source/KGeoBag.cc index 27ac11ae5..8c61ae3c9 100644 --- a/KGeoBag/Source/Application/Source/KGeoBag.cc +++ b/KGeoBag/Source/Application/Source/KGeoBag.cc @@ -15,7 +15,7 @@ int main(int argc, char** argv) } auto& tXML = KXMLInitializer::GetInstance(); - auto* tTokenizer = tXML.Configure(argc, argv, false); // process files below + auto* tTokenizer = tXML.Configure(argc, argv); // process extra files below //tXML.DumpConfiguration(); @@ -23,7 +23,7 @@ int main(int argc, char** argv) for (auto tFilename : tXML.GetArguments().ParameterList()) { - coremsg(eInfo) << "processing file <" << tFilename << "> ..." << eom; + coremsg(eNormal) << "processing file <" << tFilename << "> ..." << eom; auto* tFile = new KTextFile(); tFile->AddToNames(tFilename); tTokenizer->ProcessFile(tFile); diff --git a/KGeoBag/Source/Bindings/Visualization/Root/Include/KGROOTGeometryPainterBuilder.hh b/KGeoBag/Source/Bindings/Visualization/Root/Include/KGROOTGeometryPainterBuilder.hh index 77019b65a..e0d2edb85 100644 --- a/KGeoBag/Source/Bindings/Visualization/Root/Include/KGROOTGeometryPainterBuilder.hh +++ b/KGeoBag/Source/Bindings/Visualization/Root/Include/KGROOTGeometryPainterBuilder.hh @@ -19,6 +19,14 @@ template<> inline bool KGROOTGeometryPainterBuilder::AddAttribute(KContainer* aC aContainer->CopyTo(fObject, &KNamed::SetName); return true; } + if (aContainer->GetName() == "file") { + aContainer->CopyTo(fObject, &KGROOTGeometryPainter::SetFile); + return true; + } + if (aContainer->GetName() == "path") { + aContainer->CopyTo(fObject, &KGROOTGeometryPainter::SetPath); + return true; + } if (aContainer->GetName() == "surfaces") { if (aContainer->AsReference().size() == 0) { return true; diff --git a/KGeoBag/Source/Bindings/Visualization/Root/Source/KGROOTGeometryPainterBuilder.cc b/KGeoBag/Source/Bindings/Visualization/Root/Source/KGROOTGeometryPainterBuilder.cc index d481bc39e..a965b59ca 100644 --- a/KGeoBag/Source/Bindings/Visualization/Root/Source/KGROOTGeometryPainterBuilder.cc +++ b/KGeoBag/Source/Bindings/Visualization/Root/Source/KGROOTGeometryPainterBuilder.cc @@ -10,6 +10,8 @@ namespace katrin { STATICINT sKGROOTGeometryPainterStructure = KGROOTGeometryPainterBuilder::Attribute("name") + + KGROOTGeometryPainterBuilder::Attribute("file") + + KGROOTGeometryPainterBuilder::Attribute("path") + KGROOTGeometryPainterBuilder::Attribute("surfaces") + KGROOTGeometryPainterBuilder::Attribute("spaces") + KGROOTGeometryPainterBuilder::Attribute("plane_normal") + diff --git a/KGeoBag/Source/Core/Include/KGSpace.hh b/KGeoBag/Source/Core/Include/KGSpace.hh index 7eb05849d..1def48eca 100644 --- a/KGeoBag/Source/Core/Include/KGSpace.hh +++ b/KGeoBag/Source/Core/Include/KGSpace.hh @@ -53,6 +53,8 @@ class KGSpace : public KTagged void AddChildSpace(KGSpace* aSpace); const KGSpace* GetParent() const; + std::string GetPath() const; + const std::vector* GetBoundaries() const; const std::vector* GetChildSurfaces() const; const std::vector* GetChildSpaces() const; diff --git a/KGeoBag/Source/Core/Include/KGSurface.hh b/KGeoBag/Source/Core/Include/KGSurface.hh index 00871326d..66b428442 100644 --- a/KGeoBag/Source/Core/Include/KGSurface.hh +++ b/KGeoBag/Source/Core/Include/KGSurface.hh @@ -50,6 +50,7 @@ class KGSurface : public KTagged void Orphan(); const KGSpace* GetParent() const; + std::string GetPath() const; protected: KGSpace* fParent; diff --git a/KGeoBag/Source/Core/Source/KGSpace.cc b/KGeoBag/Source/Core/Source/KGSpace.cc index f3312a4ec..79767daec 100644 --- a/KGeoBag/Source/Core/Source/KGSpace.cc +++ b/KGeoBag/Source/Core/Source/KGSpace.cc @@ -119,6 +119,17 @@ const KGSpace* KGSpace::GetParent() const return fParent; } +std::string KGSpace::GetPath() const +{ + string tPath = GetName(); + const KGSpace* tParent = GetParent(); + while (tParent != nullptr && tParent != KGInterface::GetInstance()->Root()) { + tPath = tParent->GetName() + "/" + tPath; + tParent = tParent->GetParent(); + } + return tPath; +} + const vector* KGSpace::GetBoundaries() const { return &fBoundaries; diff --git a/KGeoBag/Source/Core/Source/KGSurface.cc b/KGeoBag/Source/Core/Source/KGSurface.cc index 3e5dc1f79..13dae4d22 100644 --- a/KGeoBag/Source/Core/Source/KGSurface.cc +++ b/KGeoBag/Source/Core/Source/KGSurface.cc @@ -69,6 +69,16 @@ const KGSpace* KGSurface::GetParent() const return fParent; } +std::string KGSurface::GetPath() const +{ + string tPath = GetName(); + const KGSpace* tParent = GetParent(); + if (tParent != nullptr && tParent != KGInterface::GetInstance()->Root()) { + tPath = tParent->GetPath() + "/" + tPath; + } + return tPath; +} + //************* //transformable //************* diff --git a/KGeoBag/Source/Visualization/Root/Include/KGROOTGeometryPainter.hh b/KGeoBag/Source/Visualization/Root/Include/KGROOTGeometryPainter.hh index 53726fd52..fcd9891cd 100644 --- a/KGeoBag/Source/Visualization/Root/Include/KGROOTGeometryPainter.hh +++ b/KGeoBag/Source/Visualization/Root/Include/KGROOTGeometryPainter.hh @@ -111,6 +111,9 @@ class KGROOTGeometryPainter : void Display() override; void Write() override; + void WriteJSON(); + void WriteSVG(); + void AddSurface(KGSurface* aSurface); void AddSpace(KGSpace* aSpace); @@ -134,7 +137,10 @@ class KGROOTGeometryPainter : K_SET(KThreeVector, PlanePoint); K_SET(bool, SwapAxis); K_GET(KThreeVector, PlaneVectorA); - K_GET(KThreeVector, PlaneVectorB) + K_GET(KThreeVector, PlaneVectorB); + + K_SET_GET(std::string, File); + K_SET_GET(std::string, Path); public: std::string GetXAxisLabel() override; diff --git a/KGeoBag/Source/Visualization/Root/Source/KGROOTGeometryPainter.cc b/KGeoBag/Source/Visualization/Root/Source/KGROOTGeometryPainter.cc index bf7e7c700..6cd27a60b 100644 --- a/KGeoBag/Source/Visualization/Root/Source/KGROOTGeometryPainter.cc +++ b/KGeoBag/Source/Visualization/Root/Source/KGROOTGeometryPainter.cc @@ -8,6 +8,7 @@ using katrin::KFile; #include #include +#include #include #include @@ -23,6 +24,8 @@ KGROOTGeometryPainter::KGROOTGeometryPainter() : fSwapAxis(false), fPlaneVectorA(0.0, 0.0, 1.0), fPlaneVectorB(1.0, 0.0, 0.0), + fFile(""), + fPath(""), fEpsilon(1.0e-10), fROOTSpaces(), fROOTSurfaces(), @@ -225,29 +228,160 @@ std::string KGROOTGeometryPainter::GetAxisLabel(KThreeVector anAxis) void KGROOTGeometryPainter::Display() { + if (fDisplayEnabled == true) { + vismsg(eInfo) << "ROOT geometry painter drawing " << fROOTSpaces.size() << " spaces" << eom; + + for (size_t i = 0; i < fROOTSpaces.size(); i++) { + fROOTSpaces.at(i)->SetLineStyle(kSolid); + fROOTSpaces.at(i)->SetLineColor(kGreen + 2); + fROOTSpaces.at(i)->SetLineWidth(1); + fROOTSpaces.at(i)->SetFillColorAlpha(kGreen + 2, 0.8); + fROOTSpaces.at(i)->Draw("F"); + } - vismsg(eInfo) << "ROOT geometry painter drawing " << fROOTSpaces.size() << " spaces" << eom; + vismsg(eInfo) << "ROOT geometry painter drawing " << fROOTSurfaces.size() << " surfaces" << eom; - for (size_t i = 0; i < fROOTSpaces.size(); i++) { - fROOTSpaces.at(i)->SetFillColor(kGreen + 2); - fROOTSpaces.at(i)->Draw("f"); + for (size_t i = 0; i < fROOTSurfaces.size(); i++) { + fROOTSurfaces.at(i)->SetLineStyle(kSolid); + fROOTSurfaces.at(i)->SetLineColor(kBlack); + fROOTSurfaces.at(i)->SetLineWidth(1); + fROOTSurfaces.at(i)->Draw(); + } } - vismsg(eInfo) << "ROOT geometry painter drawing " << fROOTSurfaces.size() << " surfaces" << eom; + Write(); // FIXME - for (size_t i = 0; i < fROOTSurfaces.size(); i++) { - fROOTSurfaces.at(i)->SetLineColor(kBlack); - fROOTSurfaces.at(i)->SetLineWidth(1); - fROOTSurfaces.at(i)->Draw(); - } return; } void KGROOTGeometryPainter::Write() { - //root write + std::cout << "WRITE: " << fWriteEnabled << std::endl; + + //if (fWriteEnabled == true) { + + WriteJSON(); + WriteSVG(); + + //} return; } +void KGROOTGeometryPainter::WriteJSON() +{ + string tFileName; + + if (fFile.length() > 0) { + if (!fPath.empty()) { + tFileName = string(fPath) + string("/") + fFile + string(".json"); + } + else { + tFileName = string(OUTPUT_DEFAULT_DIR) + string("/") + fFile + string(".json"); + } + } + else { + if (!fPath.empty()) { + tFileName = string(fPath) + string("/") + GetName() + string(".json"); + } + else { + tFileName = string(OUTPUT_DEFAULT_DIR) + string("/") + GetName() + string(".json"); + } + } + + vismsg(eInfo) << "ROOT geometry painter writing to file <" << tFileName << ">" << eom; + + ofstream json(tFileName); + + json << "{" << endl; + + json << " \"spaces\": [" << endl; + for (auto& it : fROOTSpaces) { + + json << " ["; + for (int i = 0; i < it->Size(); i++) { + json << (i > 0 ? ", " : " ") << "[" << it->GetX()[i] << "," << it->GetY()[i] << "]"; + } + json << " ]" << (it != fROOTSpaces.back() ? "," : "") << endl; + } + json << " ]," << endl; + + json << " \"surfaces\": [" << endl; + for (auto& it : fROOTSurfaces) { + json << " ["; + for (int i = 0; i < it->Size(); i++) { + json << (i > 0 ? ", " : " ") << "[" << it->GetX()[i] << "," << it->GetY()[i] << "]"; + } + json << " ]" << (it != fROOTSurfaces.back() ? "," : "") << endl; + } + json << " ]" << endl; + + json << "}" << endl; + + json.close(); +} + +void KGROOTGeometryPainter::WriteSVG() +{ + string tFileName; + + if (fFile.length() > 0) { + if (!fPath.empty()) { + tFileName = string(fPath) + string("/") + fFile + string(".svg"); + } + else { + tFileName = string(OUTPUT_DEFAULT_DIR) + string("/") + fFile + string(".svg"); + } + } + else { + if (!fPath.empty()) { + tFileName = string(fPath) + string("/") + GetName() + string(".svg"); + } + else { + tFileName = string(OUTPUT_DEFAULT_DIR) + string("/") + GetName() + string(".svg"); + } + } + + vismsg(eInfo) << "ROOT geometry painter writing to file <" << tFileName << ">" << eom; + + ofstream svg(tFileName); + + svg << "" << endl; + + double x = GetXMin(); + double y = GetYMin(); + double w = GetXMax() - GetXMin(); + double h = GetYMax() - GetYMin(); + + const string spaceColor = "green"; + const string surfaceColor = "blue"; + const double lineWidth = 0.005; + + svg << "" + << endl; + + svg << " " << endl; + for (auto& it : fROOTSpaces) { + svg << " Size(); i++) { + svg << (i > 0 ? " " : "") << it->GetX()[i] << "," << it->GetY()[i]; + } + svg << "\" style=\"fill:" << spaceColor << "; stroke:" << spaceColor << "; stroke-width:" << lineWidth << "\"/>" + << endl; + } + + svg << " " << endl; + for (auto& it : fROOTSurfaces) { + svg << " Size(); i++) { + svg << (i > 0 ? " " : "") << it->GetX()[i] << "," << it->GetY()[i]; + } + svg << "\" style=\"fill:none; stroke:" << surfaceColor << "; stroke-width:" << lineWidth << "\"/>" << endl; + } + + svg << "" << endl; + + svg.close(); +} + void KGROOTGeometryPainter::AddSurface(KGSurface* aSurface) { fSurfaces.push_back(aSurface); diff --git a/Kassiopeia/Applications/FieldCalculation/CMakeLists.txt b/Kassiopeia/Applications/FieldCalculation/CMakeLists.txt index d18693d46..e8558d72f 100644 --- a/Kassiopeia/Applications/FieldCalculation/CMakeLists.txt +++ b/Kassiopeia/Applications/FieldCalculation/CMakeLists.txt @@ -1,11 +1,13 @@ set( APPLICATIONS_FIELDS_BASENAMES SimpleElectricFieldCalculator SimpleElectricFieldCalculatorAlongZaxis + SimpleElectricFieldCalculatorOverXYplane SimpleElectricFieldCalculatorFromFile SimpleElectricFieldCalculatorSpeedTest SimpleMagneticFieldCalculator SimpleMagneticFieldCalculatorAlongZaxis + SimpleMagneticFieldCalculatorOverXYplane SimpleMagneticFieldCalculatorFromFile SimpleMagneticFieldCalculatorSpeedTest ) diff --git a/Kassiopeia/Applications/FieldCalculation/Source/SimpleElectricFieldCalculatorAlongZaxis.cxx b/Kassiopeia/Applications/FieldCalculation/Source/SimpleElectricFieldCalculatorAlongZaxis.cxx index cdc1e8d6e..b275c614c 100644 --- a/Kassiopeia/Applications/FieldCalculation/Source/SimpleElectricFieldCalculatorAlongZaxis.cxx +++ b/Kassiopeia/Applications/FieldCalculation/Source/SimpleElectricFieldCalculatorAlongZaxis.cxx @@ -76,7 +76,8 @@ int main(int argc, char** argv) continue; } - mainmsg(eNormal) << "Electric Field at position " << tPosition << " is " << tElectricField << eom; + mainmsg(eNormal) << "Electric Field at position " << tPosition << " is " << tElectricField + << " and potential is " << tPotential << eom; outFile << std::fixed << std::setprecision(16) << tPosition.X() << "\t" << tPosition.Y() << "\t" << tPosition.Z() << "\t" << tPotential << "\t" << tElectricField.X() << "\t" << tElectricField.Y() diff --git a/Kassiopeia/Applications/FieldCalculation/Source/SimpleElectricFieldCalculatorOverXYplane.cxx b/Kassiopeia/Applications/FieldCalculation/Source/SimpleElectricFieldCalculatorOverXYplane.cxx new file mode 100644 index 000000000..e2edc001e --- /dev/null +++ b/Kassiopeia/Applications/FieldCalculation/Source/SimpleElectricFieldCalculatorOverXYplane.cxx @@ -0,0 +1,99 @@ +#include "KMessage.h" +#include "KSFieldFinder.h" +#include "KSMainMessage.h" +#include "KSRootElectricField.h" +#include "KTextFile.h" +#include "KThreeVector.hh" +#include "KXMLInitializer.hh" +#include "KXMLTokenizer.hh" + + +using namespace Kassiopeia; +using namespace katrin; +using namespace KGeoBag; +using namespace std; + +int main(int argc, char** argv) +{ + if (argc < 6) { + cout + << "usage: ./SimpleElectricFieldCalculatorOverXYplane [ <...>] " + << endl; + // output_file can be "-" (-> write to terminal) + exit(-1); + } + + mainmsg(eNormal) << "starting initialization..." << eom; + + auto& tXML = KXMLInitializer::GetInstance(); + tXML.AddDefaultIncludePath(CONFIG_DEFAULT_DIR); + tXML.Configure(argc, argv, true); + + deque tParameters = tXML.GetArguments().ParameterList(); + tParameters.pop_front(); // strip off config file name + + double Z = stod(tParameters[0]); + double Rmax = stod(tParameters[1]); + double dR = stod(tParameters[2]); + + string outFileName(tParameters[3]); + ofstream outFileStream; + streambuf* outFileBuf; + if (outFileName == "-") { + outFileBuf = cout.rdbuf(); + } + else { + outFileStream.open(outFileName); + outFileBuf = outFileStream.rdbuf(); + } + ostream outFile(outFileBuf); + + mainmsg(eNormal) << "...initialization finished" << eom; + + // initialize Electric field + KSRootElectricField tRootElectricField; + + for (size_t tIndex = 4; tIndex < tParameters.size(); tIndex++) { + KSElectricField* tElectricFieldObject = getElectricField(tParameters[tIndex]); + tElectricFieldObject->Initialize(); + tRootElectricField.AddElectricField(tElectricFieldObject); + } + + KThreeVector tElectricField; + double tPotential; + + for (double x = -Rmax; x <= Rmax; x += dR) { + for (double y = -Rmax; y <= Rmax; y += dR) { + + KThreeVector tPosition(x, y, Z); + + try { + //tRootElectricField.CalculateFieldAndPotential( tPosition, 0., tField, tPotential ); // FIXME: do not use this (issue #144) + tRootElectricField.CalculateField(tPosition, 0.0, tElectricField); + tRootElectricField.CalculatePotential(tPosition, 0.0, tPotential); + } + catch (...) { + //cout << endl; + mainmsg(eWarning) << "error - cannot calculate field at position <" << tPosition << ">" << eom; + continue; + } + + mainmsg(eNormal) << "Electric Field at position " << tPosition << " is " << tElectricField + << " and potential is " << tPotential << eom; + + outFile << std::fixed << std::setprecision(16) << tPosition.X() << "\t" << tPosition.Y() << "\t" + << tPosition.Z() << "\t" << tPotential << "\t" << tElectricField.X() << "\t" << tElectricField.Y() + << "\t" << tElectricField.Z() << endl; + } + } + + for (size_t tIndex = 4; tIndex < tParameters.size(); tIndex++) { + KSElectricField* tElectricFieldObject = getElectricField(tParameters[tIndex]); + tElectricFieldObject->Deinitialize(); + } + + if (outFileStream.is_open()) + outFileStream.close(); + + return 0; +} diff --git a/Kassiopeia/Applications/FieldCalculation/Source/SimpleMagneticFieldCalculatorOverXYplane.cxx b/Kassiopeia/Applications/FieldCalculation/Source/SimpleMagneticFieldCalculatorOverXYplane.cxx new file mode 100644 index 000000000..510fcefbd --- /dev/null +++ b/Kassiopeia/Applications/FieldCalculation/Source/SimpleMagneticFieldCalculatorOverXYplane.cxx @@ -0,0 +1,94 @@ +#include "KMessage.h" +#include "KSFieldFinder.h" +#include "KSMainMessage.h" +#include "KSRootMagneticField.h" +#include "KTextFile.h" +#include "KThreeVector.hh" +#include "KXMLInitializer.hh" +#include "KXMLTokenizer.hh" + + +using namespace Kassiopeia; +using namespace katrin; +using namespace KGeoBag; +using namespace std; + +int main(int argc, char** argv) +{ + if (argc < 6) { + cout + << "usage: ./SimpleMagneticFieldCalculatorOverXYplane [ <...>] " + << endl; + // output_file can be "-" (-> write to terminal) + exit(-1); + } + + mainmsg(eNormal) << "starting initialization..." << eom; + + auto& tXML = KXMLInitializer::GetInstance(); + tXML.AddDefaultIncludePath(CONFIG_DEFAULT_DIR); + tXML.Configure(argc, argv, true); + + deque tParameters = tXML.GetArguments().ParameterList(); + tParameters.pop_front(); // strip off config file name + + double Z = stod(tParameters[0]); + double Rmax = stod(tParameters[1]); + double dR = stod(tParameters[2]); + + string outFileName(tParameters[3]); + ofstream outFileStream; + streambuf* outFileBuf; + if (outFileName == "-") { + outFileBuf = cout.rdbuf(); + } + else { + outFileStream.open(outFileName); + outFileBuf = outFileStream.rdbuf(); + } + ostream outFile(outFileBuf); + + mainmsg(eNormal) << "...initialization finished" << eom; + + // initialize magnetic field + KSRootMagneticField tRootMagneticField; + + for (size_t tIndex = 4; tIndex < tParameters.size(); tIndex++) { + KSMagneticField* tMagneticFieldObject = getMagneticField(tParameters[tIndex]); + tMagneticFieldObject->Initialize(); + tRootMagneticField.AddMagneticField(tMagneticFieldObject); + } + + KThreeVector tMagneticField; + + for (double x = -Rmax; x <= Rmax; x += dR) { + for (double y = -Rmax; y <= Rmax; y += dR) { + + KThreeVector tPosition(x, y, Z); + + try { + tRootMagneticField.CalculateField(tPosition, 0.0, tMagneticField); + } + catch (...) { + mainmsg(eWarning) << "error - cannot calculate field at position <" << tPosition << ">" << eom; + continue; + } + + mainmsg(eNormal) << "Magnetic Field at position " << tPosition << " is " << tMagneticField << eom; + + outFile << std::fixed << std::setprecision(16) << tPosition.X() << "\t" << tPosition.Y() << "\t" + << tPosition.Z() << "\t" << tMagneticField.X() << "\t" << tMagneticField.Y() << "\t" + << tMagneticField.Z() << endl; + } + } + + for (size_t tIndex = 4; tIndex < tParameters.size(); tIndex++) { + KSMagneticField* tMagneticFieldObject = getMagneticField(tParameters[tIndex]); + tMagneticFieldObject->Deinitialize(); + } + + if (outFileStream.is_open()) + outFileStream.close(); + + return 0; +} diff --git a/Kassiopeia/Applications/Simulation/CMakeLists.txt b/Kassiopeia/Applications/Simulation/CMakeLists.txt index 000577375..7ab1237d0 100644 --- a/Kassiopeia/Applications/Simulation/CMakeLists.txt +++ b/Kassiopeia/Applications/Simulation/CMakeLists.txt @@ -1,3 +1,8 @@ +set( MAIN_BASENAMES + Kassiopeia + ParticleGenerator +) + set( MAIN_LIBRARIES ${Kommon_LIBRARIES} ${Kommon_Vtk_LIBRARIES} @@ -18,10 +23,14 @@ set( MAIN_LIBRARIES KassiopeiaBindings ) -kasper_internal_include_directories() +set( MAIN_PATH + ${CMAKE_CURRENT_SOURCE_DIR}/Source +) -# target -add_executable( Kassiopeia MACOSX_BUNDLE ${CMAKE_CURRENT_SOURCE_DIR}/Source/Kassiopeia.cxx ) -target_link_libraries (Kassiopeia ${MAIN_LIBRARIES}) +kasper_internal_include_directories() +foreach( BASENAME ${MAIN_BASENAMES} ) + add_executable( ${BASENAME} ${MAIN_PATH}/${BASENAME}.cxx ) + target_link_libraries( ${BASENAME} ${MAIN_LIBRARIES} ) +endforeach( BASENAME ) -kasper_install_executables( Kassiopeia ) +kasper_install_executables( ${MAIN_BASENAMES} ) diff --git a/Kassiopeia/Applications/Simulation/Source/Kassiopeia.cxx b/Kassiopeia/Applications/Simulation/Source/Kassiopeia.cxx index 659b2ef69..1d0b8948a 100644 --- a/Kassiopeia/Applications/Simulation/Source/Kassiopeia.cxx +++ b/Kassiopeia/Applications/Simulation/Source/Kassiopeia.cxx @@ -18,13 +18,13 @@ int main(int argc, char** argv) { if (argc == 1) { cout - << "usage: ./Kassiopeia [ <...>] [ -r variable1=value1 variable2=value ... ]" + << "usage: ./Kassiopeia [ <...>] [ -r variable1=value1 variable2=value ... ]" << endl; exit(-1); } auto& tXML = KXMLInitializer::GetInstance(); - auto* tTokenizer = tXML.Configure(argc, argv, false); // process files below + auto* tTokenizer = tXML.Configure(argc, argv, true); // process extra files below //tXML.DumpConfiguration(); @@ -33,7 +33,7 @@ int main(int argc, char** argv) KToolbox::GetInstance(); for (auto tFilename : tXML.GetArguments().ParameterList()) { - mainmsg(eInfo) << "processing file <" << tFilename << "> ..." << eom; + mainmsg(eNormal) << "processing file <" << tFilename << "> ..." << eom; auto* tFile = new KTextFile(); tFile->AddToNames(tFilename); tTokenizer->ProcessFile(tFile); diff --git a/Kassiopeia/Applications/Simulation/Source/ParticleGenerator.cxx b/Kassiopeia/Applications/Simulation/Source/ParticleGenerator.cxx new file mode 100644 index 000000000..c271bb44d --- /dev/null +++ b/Kassiopeia/Applications/Simulation/Source/ParticleGenerator.cxx @@ -0,0 +1,149 @@ +#include "KMessage.h" +#include "KSMainMessage.h" +#include "KSParticleFactory.h" +#include "KSRootElectricField.h" +#include "KSRootGenerator.h" +#include "KSRootMagneticField.h" +#include "KTextFile.h" +#include "KThreeVector.hh" +#include "KToolbox.h" +#include "KXMLInitializer.hh" +#include "KXMLTokenizer.hh" + +#ifdef KEMFIELD_USE_MPI +#include "KMPIInterface.hh" +#endif + + +using namespace Kassiopeia; +using namespace katrin; +using namespace KGeoBag; +using namespace std; + +int main(int argc, char** argv) +{ + if (argc < 2) { + cout + << "usage: ./ParticleGenerator [ <...>] " + << endl; + exit(-1); + } + + KMessageTable::GetInstance().SetPrecision(16); + + mainmsg(eNormal) << "starting initialization..." << eom; + + auto& tXML = KXMLInitializer::GetInstance(); + tXML.AddDefaultIncludePath(CONFIG_DEFAULT_DIR); + tXML.Configure(argc, argv, true); + + deque tParameters = tXML.GetArguments().ParameterList(); + tParameters.pop_front(); // strip off config file name + + string outFileName(tParameters[0]); + ofstream outFileStream; + streambuf* outFileBuf; + if (outFileName == "-") { + outFileBuf = cout.rdbuf(); + } + else { + outFileStream.open(outFileName); + outFileBuf = outFileStream.rdbuf(); + } + ostream outFile(outFileBuf); + + unsigned long nEvents = stoi(tParameters[1]); + + mainmsg(eNormal) << "...initialization finished" << eom; + + // initialize generator + KSRootGenerator tRootGenerator; + + KSRootMagneticField tMagneticField; + KSRootElectricField tElectricField; + + KSParticleFactory::GetInstance().SetMagneticField(&tMagneticField); + KSParticleFactory::GetInstance().SetElectricField(&tElectricField); + + /// TODO: cmdline option for B/E fields + + outFile << std::left << std::setfill(' ') << std::setw(20) << "# index_number\t" + << "pid\t" + << "mass\t" + << "charge\t" + << "total_spin\t" + << "gyromagnetic_ratio\t" + << "time\t" + << "position_x\t" + << "position_y\t" + << "position_z\t" + << "momentum_x\t" + << "momentum_y\t" + << "momentum_z\t" + << "kinetic_energy_ev\t" + << "magnetic_field\t" + << "electric_field\t" + << "electric_potential\t" << endl; + + for (size_t tIndex = 2; tIndex < tParameters.size(); tIndex++) { + KSGenerator* tGeneratorObject = KToolbox::GetInstance().Get(tParameters[tIndex]); + + if (!tGeneratorObject) { + mainmsg(eError) << "Generator <" << tParameters[tIndex] << "> does not exist in toolbox" << eom; + continue; + } + + tGeneratorObject->Initialize(); + tRootGenerator.SetGenerator(tGeneratorObject); + + auto tEvent = new KSEvent(); + ; + tRootGenerator.SetEvent(tEvent); + + for (unsigned i = 0; i < nEvents; i++) { + + try { + tRootGenerator.ExecuteGeneration(); + } + catch (...) { + cout << endl; + mainmsg(eWarning) << "error - cannot execute generator <" << tGeneratorObject->GetName() << ">" << eom; + return 1; + } + } + + auto tParticleQueue = tEvent->GetParticleQueue(); + + mainmsg(eNormal) << "Generator <" << tGeneratorObject->GetName() << "> created <" << tParticleQueue.size() + << "> events" << eom; + + for (auto& tParticle : tParticleQueue) { + tParticle->Print(); + + if (!tParticle->IsValid()) + continue; + + outFile << std::setw(20) << std::fixed << std::setprecision(9) << tParticle->GetIndexNumber() << "\t" + << tParticle->GetPID() << "\t" << tParticle->GetMass() << "\t" << tParticle->GetCharge() << "\t" + << tParticle->GetSpinMagnitude() << "\t" << tParticle->GetGyromagneticRatio() << "\t" + << tParticle->GetTime() << "\t" << tParticle->GetX() << "\t" << tParticle->GetY() << "\t" + << tParticle->GetZ() << "\t" << tParticle->GetPX() << "\t" << tParticle->GetPY() << "\t" + << tParticle->GetPZ() << "\t" << tParticle->GetKineticEnergy_eV() << "\t" + << tParticle->GetMagneticField().Magnitude() << "\t" << tParticle->GetElectricField().Magnitude() + << "\t" << tParticle->GetElectricPotential() << "\t" << endl; + } + + mainmsg(eNormal) << eom; + + delete tEvent; + tRootGenerator.SetEvent(nullptr); + + tRootGenerator.ClearGenerator(tGeneratorObject); + tGeneratorObject->Deinitialize(); + } + + if (outFileStream.is_open()) + outFileStream.close(); + + return 0; +} diff --git a/Kassiopeia/Bindings/CMakeLists.txt b/Kassiopeia/Bindings/CMakeLists.txt index 03ea4b99c..74d83e489 100644 --- a/Kassiopeia/Bindings/CMakeLists.txt +++ b/Kassiopeia/Bindings/CMakeLists.txt @@ -203,6 +203,7 @@ if( Kassiopeia_USE_ROOT ) Visualization/Include/KSROOTTrackPainterBuilder.h Visualization/Include/KSROOTPotentialPainterBuilder.h Visualization/Include/KSROOTMagFieldPainterBuilder.h + Visualization/Include/KSROOTZonalHarmonicsPainterBuilder.h ) if( Kassiopeia_USE_VTK ) list( APPEND BINDINGS_HEADER_FILES @@ -424,6 +425,7 @@ if( Kassiopeia_USE_ROOT ) Visualization/Source/KSROOTTrackPainterBuilder.cxx Visualization/Source/KSROOTPotentialPainterBuilder.cxx Visualization/Source/KSROOTMagFieldPainterBuilder.cxx + Visualization/Source/KSROOTZonalHarmonicsPainterBuilder.cxx ) if( Kassiopeia_USE_VTK ) list( APPEND BINDINGS_SOURCE_FILES diff --git a/Kassiopeia/Bindings/Simulation/Include/KSRootBuilder.h b/Kassiopeia/Bindings/Simulation/Include/KSRootBuilder.h index bc9f1fcd0..fafbf399f 100644 --- a/Kassiopeia/Bindings/Simulation/Include/KSRootBuilder.h +++ b/Kassiopeia/Bindings/Simulation/Include/KSRootBuilder.h @@ -6,6 +6,7 @@ #include "KMagneticField.hh" #include "KSElectricKEMField.h" #include "KSMagneticKEMField.h" +#include "KSMainMessage.h" #include "KSRoot.h" #include "KSSimulation.h" #include "KToolbox.h" @@ -26,8 +27,12 @@ template<> inline bool KSRootBuilder::AddElement(KContainer* aContainer) KToolbox::GetInstance().AddContainer(*aContainer); return true; } - // legacy support for old field bindings in the tag + + /// NOTE: deprecated legacy support for old field bindings in the tag if (aContainer->Is()) { + mainmsg(eWarning) << "legacy binding for electric field <" << aContainer->GetName() + << "> is DEPRECATED - please move objects to tag" << eom; + auto* tField = new KSElectricKEMField(); tField->SetName(aContainer->GetName()); aContainer->ReleaseTo(tField, &KSElectricKEMField::SetElectricField); @@ -35,6 +40,9 @@ template<> inline bool KSRootBuilder::AddElement(KContainer* aContainer) return true; } if (aContainer->Is()) { + mainmsg(eWarning) << "legacy binding for magnetic field <" << aContainer->GetName() + << "> is DEPRECATED - please move objects to tag" << eom; + auto* tField = new KSMagneticKEMField(); tField->SetName(aContainer->GetName()); aContainer->ReleaseTo(tField, &KSMagneticKEMField::SetMagneticField); diff --git a/Kassiopeia/Bindings/Visualization/Include/KSROOTZonalHarmonicsPainterBuilder.h b/Kassiopeia/Bindings/Visualization/Include/KSROOTZonalHarmonicsPainterBuilder.h new file mode 100644 index 000000000..b46819170 --- /dev/null +++ b/Kassiopeia/Bindings/Visualization/Include/KSROOTZonalHarmonicsPainterBuilder.h @@ -0,0 +1,109 @@ +#ifndef KSROOTZONALHARMONICSPAINTERBUILDER_H +#define KSROOTZONALHARMONICSPAINTERBUILDER_H + +#include "KComplexElement.hh" +#include "KSROOTZonalHarmonicsPainter.h" +#include "KSVisualizationMessage.h" + + +using namespace Kassiopeia; +namespace katrin +{ +typedef KComplexElement KSROOTZonalHarmonicsPainterBuilder; + +template<> inline bool KSROOTZonalHarmonicsPainterBuilder::AddAttribute(KContainer* aContainer) +{ + if (aContainer->GetName() == "name") { + aContainer->CopyTo(fObject, &KNamed::SetName); + return true; + } + if (aContainer->GetName() == "x_axis") { + aContainer->CopyTo(fObject, &KSROOTZonalHarmonicsPainter::SetXAxis); + return true; + } + if (aContainer->GetName() == "y_axis") { + aContainer->CopyTo(fObject, &KSROOTZonalHarmonicsPainter::SetYAxis); + return true; + } + if (aContainer->GetName() == "electric_field") { + aContainer->CopyTo(fObject, &KSROOTZonalHarmonicsPainter::SetElectricFieldName); + return true; + } + if (aContainer->GetName() == "magnetic_field") { + aContainer->CopyTo(fObject, &KSROOTZonalHarmonicsPainter::SetMagneticFieldName); + return true; + } + if (aContainer->GetName() == "r_min") { + aContainer->CopyTo(fObject, &KSROOTZonalHarmonicsPainter::SetRmin); + return true; + } + if (aContainer->GetName() == "r_max") { + aContainer->CopyTo(fObject, &KSROOTZonalHarmonicsPainter::SetRmax); + return true; + } + if (aContainer->GetName() == "z_min") { + aContainer->CopyTo(fObject, &KSROOTZonalHarmonicsPainter::SetZmin); + return true; + } + if (aContainer->GetName() == "z_max") { + aContainer->CopyTo(fObject, &KSROOTZonalHarmonicsPainter::SetZmax); + return true; + } + if (aContainer->GetName() == "z_dist") { + aContainer->CopyTo(fObject, &KSROOTZonalHarmonicsPainter::SetZdist); + return true; + } + if (aContainer->GetName() == "r_dist") { + aContainer->CopyTo(fObject, &KSROOTZonalHarmonicsPainter::SetRdist); + return true; + } + if (aContainer->GetName() == "r_steps") { + aContainer->CopyTo(fObject, &KSROOTZonalHarmonicsPainter::SetRMaxSteps); + return true; + } + if (aContainer->GetName() == "z_steps") { + aContainer->CopyTo(fObject, &KSROOTZonalHarmonicsPainter::SetZMaxSteps); + return true; + } + if (aContainer->GetName() == "path") { + aContainer->CopyTo(fObject, &KSROOTZonalHarmonicsPainter::SetPath); + return true; + } + if (aContainer->GetName() == "file") { + aContainer->CopyTo(fObject, &KSROOTZonalHarmonicsPainter::SetFile); + return true; + } + if (aContainer->GetName() == "write") { + aContainer->CopyTo(fObject, &KSROOTZonalHarmonicsPainter::SetWriteMode); + return true; + } +#if 0 + if( aContainer->GetName() == "geometry_type" ) + { + if( aContainer->AsReference< std::string >() == "volume" ) + { + vismsg(eDebug) << "setting painter to volume write mode" << eom; + aContainer->CopyTo( fObject, &KSROOTZonalHarmonicsPainter::SetGeometryType ); + return true; + } + if( aContainer->AsReference< std::string >() == "surface" ) + { + vismsg(eDebug) << "setting painter to surface write mode" << eom; + aContainer->CopyTo( fObject, &KSROOTZonalHarmonicsPainter::SetGeometryType ); + return true; + } + return false; + } + if( aContainer->GetName() == "radial_safety_margin" ) + { + aContainer->CopyTo( fObject, &KSROOTZonalHarmonicsPainter::SetRadialSafetyMargin ); + return true; + } +#endif + + return false; +} + +} // namespace katrin + +#endif // KSROOTZONALHARMONICSPAINTERBUILDER_H diff --git a/Kassiopeia/Bindings/Visualization/Source/KSROOTZonalHarmonicsPainterBuilder.cxx b/Kassiopeia/Bindings/Visualization/Source/KSROOTZonalHarmonicsPainterBuilder.cxx new file mode 100644 index 000000000..92cc9e115 --- /dev/null +++ b/Kassiopeia/Bindings/Visualization/Source/KSROOTZonalHarmonicsPainterBuilder.cxx @@ -0,0 +1,39 @@ +#include "KSROOTZonalHarmonicsPainterBuilder.h" + +#include "KROOTPadBuilder.h" +#include "KROOTWindowBuilder.h" + +using namespace Kassiopeia; +using namespace std; + +namespace katrin +{ + +STATICINT sKSROOTZonalHarmonicsPainterStructure = + KSROOTZonalHarmonicsPainterBuilder::Attribute("name") + + KSROOTZonalHarmonicsPainterBuilder::Attribute("x_axis") + + KSROOTZonalHarmonicsPainterBuilder::Attribute("y_axis") + + KSROOTZonalHarmonicsPainterBuilder::Attribute("electric_field") + + KSROOTZonalHarmonicsPainterBuilder::Attribute("magnetic_field") + + KSROOTZonalHarmonicsPainterBuilder::Attribute("r_min") + + KSROOTZonalHarmonicsPainterBuilder::Attribute("r_max") + + KSROOTZonalHarmonicsPainterBuilder::Attribute("z_min") + + KSROOTZonalHarmonicsPainterBuilder::Attribute("z_max") + + KSROOTZonalHarmonicsPainterBuilder::Attribute("z_dist") + + KSROOTZonalHarmonicsPainterBuilder::Attribute("r_dist") + + KSROOTZonalHarmonicsPainterBuilder::Attribute("r_steps") + + KSROOTZonalHarmonicsPainterBuilder::Attribute("z_steps") + + KSROOTZonalHarmonicsPainterBuilder::Attribute("path") + + KSROOTZonalHarmonicsPainterBuilder::Attribute("file") + + KSROOTZonalHarmonicsPainterBuilder::Attribute("write"); +//KSROOTZonalHarmonicsPainterBuilder::Attribute< string >("geometry_type" ) + +//KSROOTZonalHarmonicsPainterBuilder::Attribute< double >( "radial_safety_margin" ); + + +STATICINT sKSROOTZonalHarmonicsPainterWindow = + KROOTWindowBuilder::ComplexElement("root_zh_painter"); + +STATICINT sKSROOTZonalHarmonicsPainterPad = + KROOTPadBuilder::ComplexElement("root_zh_painter"); + +} // namespace katrin diff --git a/Kassiopeia/CMakeLists.txt b/Kassiopeia/CMakeLists.txt index dc28effed..b0ac27aa3 100644 --- a/Kassiopeia/CMakeLists.txt +++ b/Kassiopeia/CMakeLists.txt @@ -13,7 +13,7 @@ endif() # Kassiopeia version set( MODULE_VERSION_MAJOR 3 ) set( MODULE_VERSION_MINOR 7 ) -set( MODULE_VERSION_PATCH 0 ) +set( MODULE_VERSION_PATCH 3 ) set( MODULE_VERSION "${MODULE_VERSION_MAJOR}.${MODULE_VERSION_MINOR}.${MODULE_VERSION_PATCH}" ) #project( Kassiopeia VERSION ${MODULE_VERSION} ) diff --git a/Kassiopeia/Visualization/CMakeLists.txt b/Kassiopeia/Visualization/CMakeLists.txt index 1aa7618b4..868d7e43f 100644 --- a/Kassiopeia/Visualization/CMakeLists.txt +++ b/Kassiopeia/Visualization/CMakeLists.txt @@ -8,6 +8,7 @@ if( Kassiopeia_USE_ROOT ) KSROOTTrackPainter.h KSROOTPotentialPainter.h KSROOTMagFieldPainter.h + KSROOTZonalHarmonicsPainter.h ) if( Kassiopeia_USE_VTK ) set( VISUALIZATION_HEADER_BASENAMES @@ -37,6 +38,7 @@ if( Kassiopeia_USE_ROOT ) KSROOTTrackPainter.cxx KSROOTPotentialPainter.cxx KSROOTMagFieldPainter.cxx + KSROOTZonalHarmonicsPainter.cxx ) if( Kassiopeia_USE_VTK ) set( VISUALIZATION_SOURCE_BASENAMES diff --git a/Kassiopeia/Visualization/Include/KSROOTZonalHarmonicsPainter.h b/Kassiopeia/Visualization/Include/KSROOTZonalHarmonicsPainter.h new file mode 100644 index 000000000..e56af813f --- /dev/null +++ b/Kassiopeia/Visualization/Include/KSROOTZonalHarmonicsPainter.h @@ -0,0 +1,88 @@ +#ifndef KSROOTZONALHARMONICSPAINTER_H +#define KSROOTZONALHARMONICSPAINTER_H + +#include "KROOTWindow.h" +using katrin::KROOTWindow; + +#include "KROOTPainter.h" +using katrin::KROOTPainter; + +#include "KField.h" +#include "TGraph.h" +#include "TPolyMarker.h" + +namespace Kassiopeia +{ + +class KSROOTZonalHarmonicsPainter : public KROOTPainter +{ + public: + KSROOTZonalHarmonicsPainter(); + ~KSROOTZonalHarmonicsPainter() override; + + void Render() override; + void Display() override; + void Write() override; + + double GetXMin() override; + double GetXMax() override; + double GetYMin() override; + double GetYMax() override; + + std::string GetXAxisLabel() override; + std::string GetYAxisLabel() override; + + private: + ; + K_SET(std::string, XAxis); + ; + K_SET(std::string, YAxis); + ; + K_SET(double, Zmin); + ; + K_SET(double, Zmax); + ; + K_SET(double, Rmin); + ; + K_SET(double, Rmax); + ; + K_SET(double, Zdist); + ; + K_SET(double, Rdist); + ; + K_SET(unsigned, ZMaxSteps); + ; + K_SET(unsigned, RMaxSteps); + ; + K_SET(std::string, ElectricFieldName); + ; + K_SET(std::string, MagneticFieldName); + ; + K_SET(std::string, File); + ; + K_SET(std::string, Path); + ; + K_SET(bool, DrawSourcePoints); + ; + K_SET(bool, DrawSourcePointArea); + //;K_SET( std::string, GeometryType ); + //;K_SET( double, RadialSafetyMargin ); + + //std::vector > fZRPoints; + + // TODO: use TPolyLine + TGraph* fElZHConvergenceGraph; + TGraph* fElZHCentralGraph; + TGraph* fElZHRemoteGraph; + + TGraph* fMagZHConvergenceGraph; + TGraph* fMagZHCentralGraph; + TGraph* fMagZHRemoteGraph; + + TPolyMarker* fElZHPoints; + TPolyMarker* fMagZHPoints; +}; + +} // namespace Kassiopeia + +#endif // KSROOTZONALHARMONICSPAINTER_H diff --git a/Kassiopeia/Visualization/Source/KSROOTZonalHarmonicsPainter.cxx b/Kassiopeia/Visualization/Source/KSROOTZonalHarmonicsPainter.cxx new file mode 100644 index 000000000..ae0de2559 --- /dev/null +++ b/Kassiopeia/Visualization/Source/KSROOTZonalHarmonicsPainter.cxx @@ -0,0 +1,453 @@ +#include "KSROOTZonalHarmonicsPainter.h" + +#include "KElectricZHFieldSolver.hh" +#include "KGElectrostaticBoundaryField.hh" +#include "KGStaticElectromagnetField.hh" +#include "KSObject.h" +#include "KSReadFileROOT.h" +#include "KSVisualizationMessage.h" +#include "KToolbox.h" +#include "KZonalHarmonicMagnetostaticFieldSolver.hh" +#include "TMultiGraph.h" + +#include +#include +#include + +using namespace KEMField; + +namespace +{ + +void MirrorGraph(TGraph* g) +{ + // used to shade area between graphs, see: https://root.cern.ch/root/html/tutorials/graphics/graphShade.C.html + if (g->GetN() <= 0) + return; + + int n = g->GetN(); + double* x = g->GetX(); + double* y = g->GetY(); + for (int i = 0; i < n; i++) { + g->SetPoint(2 * n - i, x[i], -y[i]); + } +} + +} // namespace + +namespace Kassiopeia +{ +KSROOTZonalHarmonicsPainter::KSROOTZonalHarmonicsPainter() : + fXAxis("z"), + fYAxis("r"), + fZmin(0.0), + fZmax(0.0), + fRmin(0.0), + fRmax(5.0), + fZdist(0.01), + fRdist(0.001), + fZMaxSteps(5000), + fRMaxSteps(5000), + fElectricFieldName(""), + fMagneticFieldName(""), + fFile(""), + fPath(""), + fDrawSourcePoints(true), + fDrawSourcePointArea(true) +//fGeometryType ( "surface" ), +//fRadialSafetyMargin( 0. ), +//fZRPoints () +{} +KSROOTZonalHarmonicsPainter::~KSROOTZonalHarmonicsPainter() {} + +void KSROOTZonalHarmonicsPainter::Render() +{ + bool autoAdjustZ = (fZmin >= fZmax); + if (autoAdjustZ) { + fZmin = std::numeric_limits::max(); + fZmax = std::numeric_limits::min(); + } + + fElZHConvergenceGraph = new TGraph(); + fElZHCentralGraph = new TGraph(); + fElZHRemoteGraph = new TGraph(); + + fMagZHConvergenceGraph = new TGraph(); + fMagZHCentralGraph = new TGraph(); + fMagZHRemoteGraph = new TGraph(); + + fElZHPoints = new TPolyMarker(); + fMagZHPoints = new TPolyMarker(); + + //fZRPoints.clear(); + + KElectricZHFieldSolver* tElZHSolver = nullptr; + if (!fElectricFieldName.empty()) { + vismsg(eNormal) << "Getting electric field " << fElectricFieldName << " from the toolbox" << eom; + auto* tElField = katrin::KToolbox::GetInstance().Get(fElectricFieldName); + if (tElField == nullptr) + vismsg(eError) << "No electric Field!" << eom; + vismsg(eNormal) << "Initialize electric field (again)" << eom; + tElField->Initialize(); + + //vismsg(eNormal) << "retrieve converter from electric field" << eom; + auto tElConverter = tElField->GetConverter(); + if (tElConverter.Null()) + vismsg(eError) << "Electric Field has no Converter! " << eom; + + //vismsg(eNormal) << "retrieve ZH solver from electric field" << eom; + tElZHSolver = dynamic_cast(&(*tElField->GetFieldSolver())); + if (tElZHSolver == nullptr) + vismsg(eError) << "Electric Field has no ZHSolver!" << eom; + + if (autoAdjustZ) { + auto z1 = tElZHSolver->GetParameters()->GetCentralZ1(); + auto z2 = tElZHSolver->GetParameters()->GetCentralZ2(); + if (z1 < fZmin) + fZmin = z1; + if (z2 > fZmax) + fZmax = z2; + } + + for (auto& sp : tElZHSolver->CentralSourcePoints()) { + fElZHPoints->SetNextPoint(sp.first, 0); + } + } + + KZonalHarmonicMagnetostaticFieldSolver* tMagZHSolver = nullptr; + if (!fMagneticFieldName.empty()) { + vismsg(eNormal) << "Getting magnetic field " << fMagneticFieldName << " from the toolbox" << eom; + auto* tMagField = katrin::KToolbox::GetInstance().Get(fMagneticFieldName); + if (tMagField == nullptr) + vismsg(eError) << "No magnetic Field!" << eom; + vismsg(eNormal) << "Initialize magnetic field (again)" << eom; + tMagField->Initialize(); + + //vismsg(eNormal) << "retrieve converter from magnetic field" << eom; + //auto tMagConverter = tMagField->GetConverter(); + //if ( tMagConverter.Null() ) + // vismsg(eError) << "Magnetic Field has no Converter! " << eom; + + //vismsg(eNormal) << "retrieve ZH solver from magnetic field" << eom; + tMagZHSolver = dynamic_cast(&(*tMagField->GetFieldSolver())); + if (tMagZHSolver == nullptr) + vismsg(eError) << "Magnetic Field has no ZHSolver!" << eom; + + if (autoAdjustZ) { + auto z1 = tMagZHSolver->GetParameters()->GetCentralZ1(); + auto z2 = tMagZHSolver->GetParameters()->GetCentralZ2(); + if (z1 < fZmin) + fZmin = z1; + if (z2 > fZmax) + fZmax = z2; + } + + for (auto& sp : tMagZHSolver->CentralSourcePoints()) { + fMagZHPoints->SetNextPoint(sp.first, 0); + } + } + + double tDeltaZ = fZdist; + if (ceil(fabs(fZmax - fZmin) / tDeltaZ) > fZMaxSteps) { + tDeltaZ = fabs(fZmax - fZmin) / fZMaxSteps; + } + double tDeltaR = fRdist; + if (ceil(fabs(fRmax - fRmin) / tDeltaR) > fRMaxSteps) { + tDeltaR = fabs(fRmax - fRmin) / fRMaxSteps; + } + + KThreeVector tPosition; + + vismsg(eNormal) << "ZH painter: start calculating convergence boundary from " << fZmin << " to " << fZmax << " m" + << eom; + for (double tZ = fZmin; tZ <= fZmax; tZ += fZdist) { + //vismsg( eInfo ) << "ZH painter: Z Position: " << tZ << eom; + double tR = 0; + + if (tElZHSolver != nullptr) { + for (auto& sp : tElZHSolver->CentralSourcePoints()) { + double rho = sp.second; + double dz = fabs(tZ - sp.first); + if (dz > rho) + continue; + + double h = sqrt(rho * rho - dz * dz); + if (h > tR) + tR = h; + } + + if (tR >= 0) { + //vismsg(eDebug) << " electric field sourcepoint radius at z=" << tZ << " is r=" << tR << " m" << eom; + fElZHConvergenceGraph->SetPoint(fElZHConvergenceGraph->GetN(), tZ, tR); + } + + // scan central convergence region + for (double tRho = 0; tRho <= tR; tRho += fRdist) { + if (!tElZHSolver->UseCentralExpansion(KThreeVector(0, tRho, tZ))) { + vismsg(eDebug) << " electric field central convergence limit at z=" << tZ << " is r=" << tRho + << " m" << eom; + fElZHCentralGraph->SetPoint(fElZHCentralGraph->GetN(), tZ, tRho - fRdist); + break; + } + } + // scan remote convergence region + for (double tRho = tR; tRho >= 0; tRho -= fRdist) { + if (!tElZHSolver->UseRemoteExpansion(KThreeVector(0, tRho, tZ))) { + vismsg(eDebug) << " electric field remote convergence limit at z=" << tZ << " is r=" << tRho + << " m" << eom; + fElZHRemoteGraph->SetPoint(fElZHRemoteGraph->GetN(), tZ, tRho + fRdist); + break; + } + } + } + + if (tMagZHSolver != nullptr) { + for (auto& sp : tMagZHSolver->CentralSourcePoints()) { + double rho = sp.second; + double dz = fabs(tZ - sp.first); + if (dz > rho) + continue; + + double h = sqrt(rho * rho - dz * dz); + if (h > tR) + tR = h; + } + + if (tR >= 0) { + //vismsg(eDebug) << " magnetic field sourcepoint radius at z=" << tZ << " is r=" << tR << " m" << eom; + fMagZHConvergenceGraph->SetPoint(fMagZHConvergenceGraph->GetN(), tZ, tR); + } + + // scan central convergence region + for (double tRho = 0; tRho <= tR; tRho += fRdist) { + if (!tMagZHSolver->UseCentralExpansion(KThreeVector(0, tRho, tZ))) { + vismsg(eDebug) << " magnetic field central convergence limit at z=" << tZ << " is r=" << tRho + << " m" << eom; + fMagZHCentralGraph->SetPoint(fMagZHCentralGraph->GetN(), tZ, tRho - fRdist); + break; + } + } + // scan central remote region + for (double tRho = tR; tRho >= 0; tRho -= fRdist) { + if (!tMagZHSolver->UseRemoteExpansion(KThreeVector(0, tRho, tZ))) { + vismsg(eDebug) << " magnetic field remote convergence limit at z=" << tZ << " is r=" << tRho + << " m" << eom; + fMagZHRemoteGraph->SetPoint(fMagZHRemoteGraph->GetN(), tZ, tRho + fRdist); + break; + } + } + } + } + + if (tElZHSolver != nullptr) { + // add lower part of graphs (negative R) for filled draw, points in reverse order + MirrorGraph(fElZHConvergenceGraph); + MirrorGraph(fElZHCentralGraph); + MirrorGraph(fElZHRemoteGraph); + + fElZHConvergenceGraph->SetLineColor(kRed); + fElZHConvergenceGraph->SetLineStyle(kDotted); + + fElZHCentralGraph->SetLineColor(kRed); + fElZHCentralGraph->SetFillColorAlpha(kRed, 0.1); + + fElZHRemoteGraph->SetLineColor(kRed); + + fElZHPoints->SetMarkerColor(kRed); + fElZHPoints->SetMarkerStyle(kFullDotMedium); + } + + if (tMagZHSolver != nullptr) { + // add lower part of graph (negative R) for filled draw, points in reverse order + MirrorGraph(fMagZHConvergenceGraph); + MirrorGraph(fMagZHCentralGraph); + MirrorGraph(fMagZHRemoteGraph); + + fMagZHConvergenceGraph->SetLineColor(kBlue); + fMagZHConvergenceGraph->SetLineStyle(kDotted); + + fMagZHCentralGraph->SetLineColor(kBlue); + fMagZHCentralGraph->SetFillColorAlpha(kBlue, 0.1); + + fMagZHRemoteGraph->SetLineColor(kBlue); + + fMagZHPoints->SetMarkerColor(kBlue); + fMagZHPoints->SetMarkerStyle(kFullDotMedium); + } + + return; +} + +void KSROOTZonalHarmonicsPainter::Display() +{ + if (fDisplayEnabled == true) { + // draw order may be important, who knows what ROOT does anyway ... + + auto tGraphsFilled = new TMultiGraph(); + auto tGraphs = new TMultiGraph(); + + if (fElZHCentralGraph->GetN() > 0) + tGraphsFilled->Add(fElZHCentralGraph); + if (fMagZHCentralGraph->GetN() > 0) + tGraphsFilled->Add(fMagZHCentralGraph); + + // FIXME: this is confusing + //if (fElZHRemoteGraph->GetN() > 0) + // tGraphs->Add(fElZHRemoteGraph); + //if (fMagZHRemoteGraph->GetN() > 0) + // tGraphs->Add(fMagZHRemoteGraph); + + // TODO: make this an option + if (fDrawSourcePointArea) { + if (fElZHConvergenceGraph->GetN() > 0) + tGraphs->Add(fElZHConvergenceGraph); + if (fMagZHConvergenceGraph->GetN() > 0) + tGraphs->Add(fMagZHConvergenceGraph); + } + + tGraphsFilled->Draw("fl"); + tGraphs->Draw("l"); + + // TODO: make this an option + if (fDrawSourcePoints) { + if (fElZHPoints->GetN() > 0) + fElZHPoints->Draw(); + if (fMagZHPoints->GetN() > 0) + fMagZHPoints->Draw(); + } + } + + return; +} + +void KSROOTZonalHarmonicsPainter::Write() +{ + if (fWriteEnabled == false) + return; + +#if 0 + std::string tFile; + + if( fFile.length() > 0 ) + { + if( fPath.empty() ) + { + tFile = string( OUTPUT_DEFAULT_DIR ) + string("/") + fFile; + } + else + { + tFile = fPath + string("/") + fFile; + } + } + else + { + tFile = string( OUTPUT_DEFAULT_DIR ) + string("/") + GetName() + string( ".xml" ); + } + + vismsg( eNormal ) << "writing convergence region to file <" << tFile << ">" << eom; + std::ofstream tGeometryOutput( tFile.c_str() ); + vismsg( eNormal ) << "write mode is " << fGeometryType << eom; + + if( fGeometryType == "surface" ) + { + tGeometryOutput << "" << "\n"; + tGeometryOutput << " " << "\n"; + tGeometryOutput << " " << "\n"; + tGeometryOutput << " " << "\n"; + tGeometryOutput << " " << "\n"; + + for( unsigned int i=1; i" << "\n"; + } + + tGeometryOutput << " " << "\n"; + tGeometryOutput << " " << "\n"; + tGeometryOutput << " " << "\n"; + tGeometryOutput << "" << "\n"; + tGeometryOutput.close(); + + return; + + } else if ( fGeometryType == "volume" ) { + + tGeometryOutput << "" << "\n"; + tGeometryOutput << " " << "\n"; + + for( unsigned int i=1; i" << "\n"; + } else { + tGeometryOutput << " " << "\n"; + } + + } + + tGeometryOutput << " " << "\n"; + tGeometryOutput << "" << "\n"; + + for( unsigned int i=1; i\n"; + } + + tGeometryOutput << "\n"; + tGeometryOutput << "" << "\n"; + + tGeometryOutput.close(); + return; + + } else { + vismsg(eWarning) << "wrong fGeometryType for root_zh_painter, cannot write output!" << eom; + } + + tGeometryOutput.close(); +#endif + + return; +} + +double KSROOTZonalHarmonicsPainter::GetXMin() +{ + double tMin(std::numeric_limits::max()); + return tMin; +} +double KSROOTZonalHarmonicsPainter::GetXMax() +{ + double tMax(-1.0 * std::numeric_limits::max()); + return tMax; +} + +double KSROOTZonalHarmonicsPainter::GetYMin() +{ + double tMin(std::numeric_limits::max()); + return tMin; +} +double KSROOTZonalHarmonicsPainter::GetYMax() +{ + double tMax(-1.0 * std::numeric_limits::max()); + return tMax; +} + +std::string KSROOTZonalHarmonicsPainter::GetXAxisLabel() +{ + return fXAxis; +} + +std::string KSROOTZonalHarmonicsPainter::GetYAxisLabel() +{ + return fYAxis; +} + +} // namespace Kassiopeia diff --git a/Kassiopeia/XML/CMakeLists.txt b/Kassiopeia/XML/CMakeLists.txt index 1765deb2c..68eba68c3 100644 --- a/Kassiopeia/XML/CMakeLists.txt +++ b/Kassiopeia/XML/CMakeLists.txt @@ -42,3 +42,8 @@ else( Kassiopeia_USE_VTK ) Examples/PhotoMultiplierTubeSimulation.xml ) endif( Kassiopeia_USE_VTK ) + +kasper_install_config_subdir( Examples + Examples/Fields.xml + Examples/Generators.xml +) diff --git a/Kassiopeia/XML/Examples/AnalyticSimulation.xml b/Kassiopeia/XML/Examples/AnalyticSimulation.xml index 7448b138d..cf53a5860 100644 --- a/Kassiopeia/XML/Examples/AnalyticSimulation.xml +++ b/Kassiopeia/XML/Examples/AnalyticSimulation.xml @@ -68,18 +68,18 @@ @@ -183,7 +183,8 @@ - + @@ -222,18 +223,18 @@ diff --git a/Kassiopeia/XML/Examples/DipoleTrapMeshedSpaceSimulation.xml b/Kassiopeia/XML/Examples/DipoleTrapMeshedSpaceSimulation.xml index 8a2e8970f..386b99875 100644 --- a/Kassiopeia/XML/Examples/DipoleTrapMeshedSpaceSimulation.xml +++ b/Kassiopeia/XML/Examples/DipoleTrapMeshedSpaceSimulation.xml @@ -42,12 +42,12 @@ @@ -55,13 +55,13 @@ @@ -70,13 +70,13 @@ @@ -85,13 +85,13 @@ @@ -122,19 +122,19 @@ - - - + + + - - - + + + - - - + + + @@ -191,62 +191,62 @@ + name="field_electromagnet" + directory="[KEMFIELD_CACHE]" + file="DipoleTrapMagnets.kbd" + system="world/dipole_trap" + spaces="world/dipole_trap/@magnet_tag" + > + name="field_electrostatic" + directory="[KEMFIELD_CACHE]" + file="DipoleTrapElectrodes.kbd" + system="world/dipole_trap" + surfaces="world/dipole_trap/@electrode_tag" + symmetry="axial" + > @@ -284,7 +284,8 @@ - + @@ -347,7 +348,8 @@ - + @@ -402,43 +404,43 @@ - + - + diff --git a/Kassiopeia/XML/Examples/DipoleTrapSimulation.xml b/Kassiopeia/XML/Examples/DipoleTrapSimulation.xml index b08fc9584..1153e8729 100644 --- a/Kassiopeia/XML/Examples/DipoleTrapSimulation.xml +++ b/Kassiopeia/XML/Examples/DipoleTrapSimulation.xml @@ -40,12 +40,12 @@ @@ -53,13 +53,13 @@ @@ -67,13 +67,13 @@ @@ -146,63 +146,63 @@ + name="field_electromagnet" + directory="[KEMFIELD_CACHE]" + file="DipoleTrapMagnets.kbd" + system="world/dipole_trap" + spaces="world/dipole_trap/@magnet_tag" + > + name="field_electrostatic" + directory="[KEMFIELD_CACHE]" + file="DipoleTrapElectrodes.kbd" + system="world/dipole_trap" + surfaces="world/dipole_trap/@electrode_tag" + symmetry="axial" + > @@ -272,7 +272,8 @@ - + @@ -326,18 +327,18 @@ diff --git a/Kassiopeia/XML/Examples/Fields.xml b/Kassiopeia/XML/Examples/Fields.xml new file mode 100644 index 000000000..e9690ec3d --- /dev/null +++ b/Kassiopeia/XML/Examples/Fields.xml @@ -0,0 +1,326 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Kassiopeia/XML/Examples/Generators.xml b/Kassiopeia/XML/Examples/Generators.xml new file mode 100644 index 000000000..1b1b3f3b1 --- /dev/null +++ b/Kassiopeia/XML/Examples/Generators.xml @@ -0,0 +1,316 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Kassiopeia/XML/Examples/PhotoMultiplierTubeSimulation.xml b/Kassiopeia/XML/Examples/PhotoMultiplierTubeSimulation.xml index b740c169f..bb1efe21f 100644 --- a/Kassiopeia/XML/Examples/PhotoMultiplierTubeSimulation.xml +++ b/Kassiopeia/XML/Examples/PhotoMultiplierTubeSimulation.xml @@ -47,62 +47,89 @@ - + - + - + - - - - - - - - - - - + + + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - - + + + + + + + + + + @@ -137,33 +164,41 @@ - - - + + + - - - + + + - - - + + + - - - + + + @@ -223,8 +258,8 @@ - - + + @@ -237,24 +272,41 @@ - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + - + @@ -263,7 +315,8 @@ - + @@ -278,82 +331,82 @@ + name="kemfield_e" + directory="[KEMFIELD_CACHE]" + file="PhotoMultiplierTubeElectrodes.kbd" + system="world/photomultiplier_tube" + surfaces="world/photomultiplier_tube/@electrode_tag" + symmetry="none" + minimum_element_area="1e-15" + maximum_element_aspect_ratio="50" + > - + solver_name="gmres" + preconditioner="implicit_krylov" + tolerance="1e-6" + preconditioner_tolerance="0.01" + max_iterations="500" + max_preconditioner_iterations="30" + iterations_between_restarts="250" + preconditioner_degree="1" + intermediate_save_interval="1" + use_display="true" + show_plot="false" + use_timer="true" + > + + + - - @@ -396,10 +449,10 @@ - + required_energy_per_particle_ev="10."/> @@ -413,9 +466,9 @@ @@ -529,17 +582,17 @@ diff --git a/Kassiopeia/XML/Examples/QuadrupoleTrapSimulation.xml b/Kassiopeia/XML/Examples/QuadrupoleTrapSimulation.xml index a5cb5d25c..0f121285a 100644 --- a/Kassiopeia/XML/Examples/QuadrupoleTrapSimulation.xml +++ b/Kassiopeia/XML/Examples/QuadrupoleTrapSimulation.xml @@ -97,7 +97,7 @@ - + @@ -144,56 +144,56 @@ + name="field_electromagnet" + directory="[KEMFIELD_CACHE]" + file="QuadrupoleTrapMagnets.kbd" + system="world/quadrupole_trap" + spaces="world/quadrupole_trap/@magnet_tag" + > + name="field_electrostatic" + directory="[KEMFIELD_CACHE]" + file="QuadrupoleTrapElectrodes.kbd" + system="world/quadrupole_trap" + surfaces="world/quadrupole_trap/@electrode_tag" + symmetry="axial" + > @@ -341,7 +341,8 @@ - + @@ -385,18 +386,18 @@ diff --git a/Kassiopeia/XML/Examples/ToricTrapSimulation.xml b/Kassiopeia/XML/Examples/ToricTrapSimulation.xml index 7b789e2dd..8b5a45a31 100644 --- a/Kassiopeia/XML/Examples/ToricTrapSimulation.xml +++ b/Kassiopeia/XML/Examples/ToricTrapSimulation.xml @@ -42,12 +42,12 @@ @@ -61,10 +61,10 @@ @@ -123,15 +123,15 @@ @@ -152,38 +152,38 @@ + name="field_electromagnet" + directory="[KEMFIELD_CACHE]" + file="ToricTrapMagnets.kbd" + system="world/toric_trap" + spaces="world/toric_trap/@magnet_tag" + > @@ -283,7 +283,8 @@ - + @@ -335,18 +336,18 @@ diff --git a/Kassiopeia/XML/Examples/VTK/AnalyticSimulation.xml b/Kassiopeia/XML/Examples/VTK/AnalyticSimulation.xml index c3ecd8341..63820ebb3 100644 --- a/Kassiopeia/XML/Examples/VTK/AnalyticSimulation.xml +++ b/Kassiopeia/XML/Examples/VTK/AnalyticSimulation.xml @@ -3,7 +3,7 @@ - + @@ -35,7 +35,7 @@ - + @@ -47,13 +47,13 @@ - + - + @@ -68,18 +68,18 @@ - + @@ -100,7 +100,7 @@ - + @@ -111,9 +111,9 @@ - + - + @@ -135,7 +135,7 @@ - + @@ -149,7 +149,7 @@ - + @@ -183,7 +183,8 @@ - + @@ -222,50 +223,50 @@ - - + diff --git a/Kassiopeia/XML/Examples/VTK/DipoleTrapMeshedSpaceSimulation.xml b/Kassiopeia/XML/Examples/VTK/DipoleTrapMeshedSpaceSimulation.xml index 7a1e313bd..2990061ae 100644 --- a/Kassiopeia/XML/Examples/VTK/DipoleTrapMeshedSpaceSimulation.xml +++ b/Kassiopeia/XML/Examples/VTK/DipoleTrapMeshedSpaceSimulation.xml @@ -42,12 +42,12 @@ @@ -55,13 +55,13 @@ @@ -70,13 +70,13 @@ @@ -85,13 +85,13 @@ @@ -122,19 +122,19 @@ - - - + + + - - - + + + - - - + + + @@ -191,63 +191,63 @@ + name="field_electromagnet" + directory="[KEMFIELD_CACHE]" + file="DipoleTrapMagnets.kbd" + system="world/dipole_trap" + spaces="world/dipole_trap/@magnet_tag" + > + name="field_electrostatic" + directory="[KEMFIELD_CACHE]" + file="DipoleTrapElectrodes.kbd" + system="world/dipole_trap" + surfaces="world/dipole_trap/@electrode_tag" + symmetry="axial" + > @@ -285,7 +285,8 @@ - + @@ -349,7 +350,8 @@ - + @@ -408,76 +410,76 @@ - + - + + name="vtk_window" + enable_display="true" + enable_write="true" + frame_title="KGeoBag Visualization" + frame_size_x="1024" + frame_size_y="768" + frame_color_red=".2" + frame_color_green=".2" + frame_color_blue=".2" + view_angle="45" + eye_angle="0.5" + multi_samples="4" + depth_peeling="10" + > diff --git a/Kassiopeia/XML/Examples/VTK/DipoleTrapSimulation.xml b/Kassiopeia/XML/Examples/VTK/DipoleTrapSimulation.xml index 49876ad5d..ac157da1c 100644 --- a/Kassiopeia/XML/Examples/VTK/DipoleTrapSimulation.xml +++ b/Kassiopeia/XML/Examples/VTK/DipoleTrapSimulation.xml @@ -40,12 +40,12 @@ @@ -53,13 +53,13 @@ @@ -67,13 +67,13 @@ @@ -146,64 +146,64 @@ + name="field_electromagnet" + directory="[KEMFIELD_CACHE]" + file="DipoleTrapMagnets.kbd" + system="world/dipole_trap" + spaces="world/dipole_trap/@magnet_tag" + > + name="field_electrostatic" + directory="[KEMFIELD_CACHE]" + file="DipoleTrapElectrodes.kbd" + system="world/dipole_trap" + surfaces="world/dipole_trap/@electrode_tag" + symmetry="axial" + > @@ -275,7 +275,8 @@ - + @@ -333,66 +334,66 @@ + name="vtk_window" + enable_display="true" + enable_write="true" + frame_title="KGeoBag Visualization" + frame_size_x="1024" + frame_size_y="768" + frame_color_red=".2" + frame_color_green=".2" + frame_color_blue=".2" + view_angle="45" + eye_angle="0.5" + multi_samples="4" + depth_peeling="10" + > diff --git a/Kassiopeia/XML/Examples/VTK/PhotoMultiplierTubeSimulation.xml b/Kassiopeia/XML/Examples/VTK/PhotoMultiplierTubeSimulation.xml index 28b98e7bc..4153ced3c 100644 --- a/Kassiopeia/XML/Examples/VTK/PhotoMultiplierTubeSimulation.xml +++ b/Kassiopeia/XML/Examples/VTK/PhotoMultiplierTubeSimulation.xml @@ -47,62 +47,89 @@ - + - + - + - - - - - - - - - - - + + + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - - + + + + + + + + + + @@ -137,33 +164,41 @@ - - - + + + - - - + + + - - - + + + - - - + + + @@ -223,8 +258,8 @@ - - + + @@ -237,24 +272,41 @@ - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + - + @@ -263,7 +315,8 @@ - + @@ -278,82 +331,82 @@ + name="kemfield_e" + directory="[KEMFIELD_CACHE]" + file="PhotoMultiplierTubeElectrodes.kbd" + system="world/photomultiplier_tube" + surfaces="world/photomultiplier_tube/@electrode_tag" + symmetry="none" + minimum_element_area="1e-15" + maximum_element_aspect_ratio="50" + > - + solver_name="gmres" + preconditioner="implicit_krylov" + tolerance="1e-6" + preconditioner_tolerance="0.01" + max_iterations="500" + max_preconditioner_iterations="30" + iterations_between_restarts="250" + preconditioner_degree="1" + intermediate_save_interval="1" + use_display="true" + show_plot="false" + use_timer="true" + > + + + - - @@ -397,10 +450,10 @@ - + required_energy_per_particle_ev="30."/> @@ -413,16 +466,16 @@ @@ -539,50 +592,50 @@ + name="vtk_window" + enable_display="true" + enable_write="true" + frame_title="KGeoBag Visualization" + frame_size_x="1024" + frame_size_y="768" + frame_color_red=".2" + frame_color_green=".2" + frame_color_blue=".2" + view_angle="45" + eye_angle="0.5" + multi_samples="4" + depth_peeling="10" + > diff --git a/Kassiopeia/XML/Examples/VTK/QuadrupoleTrapSimulation.xml b/Kassiopeia/XML/Examples/VTK/QuadrupoleTrapSimulation.xml index 35c8fb389..bc30f78d5 100644 --- a/Kassiopeia/XML/Examples/VTK/QuadrupoleTrapSimulation.xml +++ b/Kassiopeia/XML/Examples/VTK/QuadrupoleTrapSimulation.xml @@ -97,7 +97,7 @@ - + @@ -144,57 +144,57 @@ + name="field_electromagnet" + directory="[KEMFIELD_CACHE]" + file="QuadrupoleTrapMagnets.kbd" + system="world/quadrupole_trap" + spaces="world/quadrupole_trap/@magnet_tag" + > + name="field_electrostatic" + directory="[KEMFIELD_CACHE]" + file="QuadrupoleTrapElectrodes.kbd" + system="world/quadrupole_trap" + surfaces="world/quadrupole_trap/@electrode_tag" + symmetry="axial" + > @@ -341,7 +341,8 @@ - + @@ -389,50 +390,50 @@ + name="vtk_window" + enable_display="true" + enable_write="true" + frame_title="KGeoBag Visualization" + frame_size_x="1024" + frame_size_y="768" + frame_color_red=".2" + frame_color_green=".2" + frame_color_blue=".2" + view_angle="45" + eye_angle="0.5" + multi_samples="4" + depth_peeling="10" + > diff --git a/Kassiopeia/XML/Examples/VTK/ToricTrapSimulation.xml b/Kassiopeia/XML/Examples/VTK/ToricTrapSimulation.xml index 2f1b31423..f3e3cb8f2 100644 --- a/Kassiopeia/XML/Examples/VTK/ToricTrapSimulation.xml +++ b/Kassiopeia/XML/Examples/VTK/ToricTrapSimulation.xml @@ -42,12 +42,12 @@ @@ -61,10 +61,10 @@ @@ -123,15 +123,15 @@ @@ -148,31 +148,31 @@ + name="vtk_window" + enable_display="true" + enable_write="true" + frame_title="KGeoBag Visualization" + frame_size_x="1024" + frame_size_y="768" + frame_color_red=".2" + frame_color_green=".2" + frame_color_blue=".2" + view_angle="45" + eye_angle="0.5" + multi_samples="4" + depth_peeling="10" + > @@ -181,38 +181,38 @@ + name="field_electromagnet" + directory="[KEMFIELD_CACHE]" + file="ToricTrapMagnets.kbd" + system="world/toric_trap" + spaces="world/toric_trap/@magnet_tag" + > @@ -314,7 +314,8 @@ - + @@ -370,19 +371,19 @@ diff --git a/Kassiopeia/XML/Validation/TestRampedField.xml b/Kassiopeia/XML/Validation/TestRampedField.xml index d18bcaf67..b24a4d7df 100644 --- a/Kassiopeia/XML/Validation/TestRampedField.xml +++ b/Kassiopeia/XML/Validation/TestRampedField.xml @@ -1,39 +1,9 @@ - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + @@ -43,19 +13,19 @@ - + - - - - - - + + - + - - + - + diff --git a/Kommon/Core/Initialization/KElementBase.cc b/Kommon/Core/Initialization/KElementBase.cc index 72222d54d..795ec5e9c 100644 --- a/Kommon/Core/Initialization/KElementBase.cc +++ b/Kommon/Core/Initialization/KElementBase.cc @@ -22,19 +22,30 @@ void KElementBase::ProcessToken(KBeginElementToken* aToken) if ((fElementDepth == 0) && (fAttributeDepth == 0)) { //look up constructor method in the map, complain and exit if not found auto It = fElements->find(aToken->GetValue()); + + // show suggestions for other elements if no match was found if (It == fElements->end()) { - // show only elements with similar beginning, and limit to maximum of 10 entries + // show only elements with similar beginning, and limit to maximum of 12 entries + const unsigned tMaxElements = 12; std::string token = aToken->GetValue(); std::set elements; while (token.length() > 0) { for (auto& it : GetElements()) { - if (elements.size() >= 10) + if (elements.size() >= tMaxElements) break; if (it.substr(0, token.length()) == token) elements.insert(it); } token.resize(token.size() - 1); } + // add other elements in alphabetical order, if list still has less than 12 entries + if (elements.size() < fElements->size()) { + for (auto& it : GetElements()) { + if (elements.size() >= tMaxElements) + break; + elements.insert(it); + } + } initmsg(eWarning) << "nothing registered for element <" << aToken->GetValue() << "> in element <" << GetName() << ">" << ret; @@ -42,7 +53,7 @@ void KElementBase::ProcessToken(KBeginElementToken* aToken) << aToken->GetLine() << "> at column <" << aToken->GetColumn() << ">" << ret; initmsg(eWarning) << "available elements: " << elements; if (elements.size() < fElements->size()) - initmsg(eWarning) << " and " << (fElements->size() - elements.size()) << " more"; + initmsg(eWarning) << " + " << (fElements->size() - elements.size()) << " more (not shown)"; initmsg(eWarning) << eom; fChild = nullptr; diff --git a/Kommon/Core/Initialization/KXMLInitializer.cc b/Kommon/Core/Initialization/KXMLInitializer.cc index 5411bca75..1f22f87e4 100644 --- a/Kommon/Core/Initialization/KXMLInitializer.cc +++ b/Kommon/Core/Initialization/KXMLInitializer.cc @@ -275,8 +275,15 @@ KXMLTokenizer* KXMLInitializer::Configure(int argc, char** argv, bool processCon KMessageTable::GetInstance().SetShowShutdownMessage(); pair tConfig; - if (processConfig) + if (processConfig) { tConfig = GetConfigFile(); + } + else { + string configLocation = fArguments.GetOption("config"); + if (!configLocation.empty()) { + KINFO("Config file ignored: " << configLocation); + } + } SetupProcessChain(fArguments.OptionTable(), tConfig.first); diff --git a/Kommon/Root/Utility/KROOTWindow.cxx b/Kommon/Root/Utility/KROOTWindow.cxx index 70e296cc1..d24775203 100644 --- a/Kommon/Root/Utility/KROOTWindow.cxx +++ b/Kommon/Root/Utility/KROOTWindow.cxx @@ -162,6 +162,8 @@ void KROOTWindow::Display() void KROOTWindow::Write() { + utilmsg(eNormal) << "KROOTWindow starts to write!" << eom; + if (fWriteEnabled) { string tOutputStringRoot; string tOutputStringPNG; @@ -200,6 +202,7 @@ void KROOTWindow::Write() } } + utilmsg(eNormal) << "KROOTWindow finished to write!" << eom; return; }