From d912b95c27c8362f1df4af6e13c1a266039dfa75 Mon Sep 17 00:00:00 2001
From: Walter Van Herck <w.van.herck@fz-juelich.de>
Date: Tue, 19 Jun 2018 12:04:39 +0200
Subject: [PATCH] Move creation of averaged multilayer out of DWBAComputation
---
Core/Computation/DWBAComputation.cpp | 47 ++----------
Core/Computation/DWBAComputation.h | 6 --
Core/Computation/DWBASingleComputation.cpp | 12 +--
Core/Computation/DWBASingleComputation.h | 4 +-
Core/Computation/IComputationUtils.cpp | 89 +++++++++++++++++++++-
Core/Computation/IComputationUtils.h | 10 +++
Core/Material/MaterialFactoryFuncs.cpp | 2 +-
Core/Material/MaterialFactoryFuncs.h | 2 +-
Tests/UnitTests/Core/Other/MaterialTest.h | 8 +-
9 files changed, 114 insertions(+), 66 deletions(-)
diff --git a/Core/Computation/DWBAComputation.cpp b/Core/Computation/DWBAComputation.cpp
index d49425f498b..a8a1ecbbe26 100644
--- a/Core/Computation/DWBAComputation.cpp
+++ b/Core/Computation/DWBAComputation.cpp
@@ -75,46 +75,13 @@ void DWBAComputation::runProtected()
}
}
-std::unique_ptr<MultiLayer> DWBAComputation::getAveragedMultilayer() const
-{
- std::map<size_t, std::vector<HomogeneousRegion>> region_map;
- m_single_computation.mergeRegionMap(region_map);
- std::unique_ptr<MultiLayer> P_result(mP_multi_layer->clone());
- auto last_layer_index = P_result->numberOfLayers()-1;
- for (auto& entry : region_map)
- {
- auto i_layer = entry.first;
- if (i_layer==0 || i_layer==last_layer_index)
- continue; // skip semi-infinite layers
- auto layer_mat = P_result->layerMaterial(i_layer);
- if (!checkRegions(entry.second))
- throw std::runtime_error("DWBAComputation::getAveragedMultilayer: "
- "total volumetric fraction of particles exceeds 1!");
- auto new_mat = createAveragedMaterial(layer_mat, entry.second);
- P_result->setLayerMaterial(i_layer, new_mat);
- }
- return P_result;
-}
-
-std::unique_ptr<MultiLayer> DWBAComputation::getMultilayerForFresnel() const
-{
- std::unique_ptr<MultiLayer> P_result = m_sim_options.useAvgMaterials()
- ? getAveragedMultilayer()
- : std::unique_ptr<MultiLayer>(mP_multi_layer->clone());
- P_result->initBFields();
- return P_result;
-}
-
void DWBAComputation::initFresnelMap()
{
- auto multilayer = getMultilayerForFresnel();
- mP_fresnel_map->setMultilayer(*multilayer);
-}
-
-bool DWBAComputation::checkRegions(const std::vector<HomogeneousRegion>& regions) const
-{
- double total_fraction = 0;
- for (auto& region : regions)
- total_fraction += region.m_volume;
- return (total_fraction >= 0 && total_fraction <= 1);
+ auto region_map = m_single_computation.regionMap();
+ std::unique_ptr<MultiLayer> P_multilayer =
+ m_sim_options.useAvgMaterials()
+ ? IComputationUtils::CreateAveragedMultilayer(*mP_multi_layer, region_map)
+ : std::unique_ptr<MultiLayer>(mP_multi_layer->clone());
+ P_multilayer->initBFields();
+ mP_fresnel_map->setMultilayer(*P_multilayer);
}
diff --git a/Core/Computation/DWBAComputation.h b/Core/Computation/DWBAComputation.h
index 86b9d92e5d0..34d76386d54 100644
--- a/Core/Computation/DWBAComputation.h
+++ b/Core/Computation/DWBAComputation.h
@@ -41,14 +41,8 @@ public:
private:
void runProtected() override;
- // creates a multilayer that contains averaged materials, for use in Fresnel calculations
- std::unique_ptr<MultiLayer> getAveragedMultilayer() const;
- // creates a multilayer for use in Fresnel calculations; if needed, it calculates average
- // materials and precalculates the magnetic B fields
- std::unique_ptr<MultiLayer> getMultilayerForFresnel() const;
// sets the correct layer materials for the Fresnel map to use
void initFresnelMap();
- bool checkRegions(const std::vector<HomogeneousRegion>& regions) const;
std::vector<SimulationElement>::iterator m_begin_it, m_end_it; //!< these iterators define the span of detector bins this simulation will work on
std::unique_ptr<IFresnelMap> mP_fresnel_map; //!< Contains the information, necessary to calculate the Fresnel coefficients.
diff --git a/Core/Computation/DWBASingleComputation.cpp b/Core/Computation/DWBASingleComputation.cpp
index 4d36aeee67a..45769b32ccf 100644
--- a/Core/Computation/DWBASingleComputation.cpp
+++ b/Core/Computation/DWBASingleComputation.cpp
@@ -33,6 +33,7 @@ void DWBASingleComputation::setProgressHandler(ProgressHandler* p_progress)
void DWBASingleComputation::addLayoutComputation(ParticleLayoutComputation* p_layout_comp)
{
m_layout_comps.emplace_back(p_layout_comp);
+ p_layout_comp->mergeRegionMap(m_region_map);
}
void DWBASingleComputation::setRoughnessComputation(RoughMultiLayerComputation* p_roughness_comp)
@@ -61,14 +62,7 @@ void DWBASingleComputation::compute(SimulationElement& elem) const
}
}
-void DWBASingleComputation::mergeRegionMap(
- std::map<size_t, std::vector<HomogeneousRegion> >& region_map) const
+const std::map<size_t, std::vector<HomogeneousRegion>>& DWBASingleComputation::regionMap() const
{
- for (auto& layout_comp : m_layout_comps) {
- layout_comp->mergeRegionMap(region_map);
- }
+ return m_region_map;
}
-
-static_assert(std::is_move_constructible<DWBASingleComputation>::value,
- "DWBASingleComputation should be move constructable");
-
diff --git a/Core/Computation/DWBASingleComputation.h b/Core/Computation/DWBASingleComputation.h
index fbec6519378..4e5fc1de414 100644
--- a/Core/Computation/DWBASingleComputation.h
+++ b/Core/Computation/DWBASingleComputation.h
@@ -48,8 +48,8 @@ public:
void setSpecularBinComputation(GISASSpecularComputation* p_spec_comp);
void compute(SimulationElement& elem) const;
- //! Merges its region map into the given one (notice non-const reference parameter)
- void mergeRegionMap(std::map<size_t, std::vector<HomogeneousRegion>>& region_map) const;
+ //! Retrieves a map of regions for the calculation of averaged layers
+ const std::map<size_t, std::vector<HomogeneousRegion>>& regionMap() const;
private:
std::vector<std::unique_ptr<ParticleLayoutComputation>> m_layout_comps;
diff --git a/Core/Computation/IComputationUtils.cpp b/Core/Computation/IComputationUtils.cpp
index 9f877e78dc5..138a573e47c 100644
--- a/Core/Computation/IComputationUtils.cpp
+++ b/Core/Computation/IComputationUtils.cpp
@@ -13,19 +13,27 @@
// ************************************************************************** //
#include "IComputationUtils.h"
+#include "ILayout.h"
+#include "IParticle.h"
+#include "Layer.h"
+#include "MaterialFactoryFuncs.h"
#include "MatrixFresnelMap.h"
#include "MultiLayer.h"
#include "ScalarFresnelMap.h"
#include "SimulationOptions.h"
+#include "SlicedFormFactorList.h"
-namespace IComputationUtils {
+namespace {
+void ScaleRegionMap(std::map<size_t, std::vector<HomogeneousRegion>>& region_map, double factor);
+bool CheckRegions(const std::vector<HomogeneousRegion>& regions);
+}
+namespace IComputationUtils {
std::unique_ptr<IFresnelMap> CreateFresnelMap(const MultiLayer& multilayer,
const SimulationOptions& sim_options,
bool allow_average_layers)
{
std::unique_ptr<IFresnelMap> P_result;
-
if (!multilayer.requiresMatrixRTCoefficients())
P_result.reset(new ScalarFresnelMap());
else
@@ -37,8 +45,83 @@ std::unique_ptr<IFresnelMap> CreateFresnelMap(const MultiLayer& multilayer,
if (!allow_average_layers && sim_options.useAvgMaterials())
throw std::runtime_error("Error in IComputationUtils::CreateFresnelMap: using averaged "
"materials is not allowed for this computation");
-
P_result->setMultilayer(multilayer);
return P_result;
}
+
+std::unique_ptr<MultiLayer> CreateAveragedMultilayer(
+ const MultiLayer& multilayer,
+ const std::map<size_t, std::vector<HomogeneousRegion>>& region_map)
+{
+ auto local_region_map = region_map.size() == 0
+ ? GetRegionMap(multilayer)
+ : region_map;
+ std::unique_ptr<MultiLayer> P_result(multilayer.clone());
+ auto last_layer_index = P_result->numberOfLayers()-1;
+ for (auto& entry : local_region_map)
+ {
+ auto i_layer = entry.first;
+ if (i_layer==0 || i_layer==last_layer_index)
+ continue; // skip semi-infinite layers
+ auto layer_mat = P_result->layerMaterial(i_layer);
+ if (!CheckRegions(entry.second))
+ throw std::runtime_error("IComputationUtils::CreateAveragedMultilayer: "
+ "total volumetric fraction of particles exceeds 1!");
+ auto new_mat = CreateAveragedMaterial(layer_mat, entry.second);
+ P_result->setLayerMaterial(i_layer, new_mat);
+ }
+ return P_result;
+}
+
+std::map<size_t, std::vector<HomogeneousRegion>> GetRegionMap(const MultiLayer& multilayer)
+{
+ std::map<size_t, std::vector<HomogeneousRegion>> result_map;
+ size_t nLayers = multilayer.numberOfLayers();
+ for (size_t i=0; i<nLayers; ++i) {
+ const Layer* layer = multilayer.layer(i);
+ for (auto p_layout : layer->layouts()) {
+ double layout_abundance = p_layout->getTotalAbundance();
+ double scale_factor = p_layout->totalParticleSurfaceDensity()/layout_abundance;
+ for (const IParticle* particle: p_layout->particles()) {
+ auto sliced_ffs = SlicedFormFactorList::CreateSlicedFormFactors(
+ *particle, multilayer, i);
+ auto particle_region_map = sliced_ffs.regionMap();
+ double particle_scaling = particle->abundance()*scale_factor;
+ ScaleRegionMap(particle_region_map, particle_scaling);
+ MergeRegionMap(result_map, particle_region_map);
+ }
+ }
+ }
+ return result_map;
}
+
+void MergeRegionMap(std::map<size_t, std::vector<HomogeneousRegion>>& dest,
+ const std::map<size_t, std::vector<HomogeneousRegion>>& source)
+{
+ for (auto& entry : source) {
+ size_t layer_index = entry.first;
+ auto regions = entry.second;
+ dest[layer_index].insert(dest[layer_index].begin(), regions.begin(), regions.end());
+ }
+}
+
+} // namespace IComputationUtils
+
+namespace {
+void ScaleRegionMap(std::map<size_t, std::vector<HomogeneousRegion>>& region_map, double factor)
+{
+ for (auto& entry : region_map) {
+ for (auto& region : entry.second) {
+ region.m_volume *= factor;
+ }
+ }
+}
+bool CheckRegions(const std::vector<HomogeneousRegion>& regions)
+{
+ double total_fraction = 0.0;
+ for (auto& region : regions)
+ total_fraction += region.m_volume;
+ return (total_fraction >= 0 && total_fraction <= 1);
+}
+} // unnamed namespace
+
diff --git a/Core/Computation/IComputationUtils.h b/Core/Computation/IComputationUtils.h
index a577286e3f0..2910068fcc6 100644
--- a/Core/Computation/IComputationUtils.h
+++ b/Core/Computation/IComputationUtils.h
@@ -15,8 +15,11 @@
#ifndef ICOMPUTATIONUTILS_H
#define ICOMPUTATIONUTILS_H
+#include <map>
#include <memory>
+#include <vector>
+class HomogeneousRegion;
class IFresnelMap;
class MultiLayer;
class SimulationOptions;
@@ -25,6 +28,13 @@ namespace IComputationUtils {
std::unique_ptr<IFresnelMap> CreateFresnelMap(const MultiLayer& multilayer,
const SimulationOptions& sim_options,
bool allow_average_layers=true);
+// creates a multilayer that contains averaged materials, for use in Fresnel calculations
+std::unique_ptr<MultiLayer> CreateAveragedMultilayer(
+ const MultiLayer& multilayer,
+ const std::map<size_t, std::vector<HomogeneousRegion>>& region_map = {});
+std::map<size_t, std::vector<HomogeneousRegion>> GetRegionMap(const MultiLayer& multilayer);
+void MergeRegionMap(std::map<size_t, std::vector<HomogeneousRegion>>& dest,
+ const std::map<size_t, std::vector<HomogeneousRegion>>& source);
}
#endif // ICOMPUTATIONUTILS_H
diff --git a/Core/Material/MaterialFactoryFuncs.cpp b/Core/Material/MaterialFactoryFuncs.cpp
index 065aad2f102..3379283153b 100644
--- a/Core/Material/MaterialFactoryFuncs.cpp
+++ b/Core/Material/MaterialFactoryFuncs.cpp
@@ -62,7 +62,7 @@ Material MaterialBySLD(const std::string& name, double sld_real, double sld_imag
return Material(std::move(mat_impl));
}
-Material createAveragedMaterial(const Material& layer_mat,
+Material CreateAveragedMaterial(const Material& layer_mat,
const std::vector<HomogeneousRegion>& regions)
{
// determine the type of returned material
diff --git a/Core/Material/MaterialFactoryFuncs.h b/Core/Material/MaterialFactoryFuncs.h
index c5f4b16f094..905984871e0 100644
--- a/Core/Material/MaterialFactoryFuncs.h
+++ b/Core/Material/MaterialFactoryFuncs.h
@@ -68,7 +68,7 @@ BA_CORE_API_ Material MaterialBySLD(const std::string& name, double sld_real, do
//! Creates averaged material. Square refractive index of returned material is arithmetic mean over
//! _regions_ and _layer_mat_. Magnetization (if present) is averaged linearly.
-BA_CORE_API_ Material createAveragedMaterial(const Material& layer_mat,
+BA_CORE_API_ Material CreateAveragedMaterial(const Material& layer_mat,
const std::vector<HomogeneousRegion>& regions);
#endif //SWIG
diff --git a/Tests/UnitTests/Core/Other/MaterialTest.h b/Tests/UnitTests/Core/Other/MaterialTest.h
index 5c2a787c15f..4987c7e3392 100644
--- a/Tests/UnitTests/Core/Other/MaterialTest.h
+++ b/Tests/UnitTests/Core/Other/MaterialTest.h
@@ -163,7 +163,7 @@ TEST_F(MaterialTest, AveragedMaterialTest)
const std::vector<HomogeneousRegion> regions
= {HomogeneousRegion{0.25, material}, HomogeneousRegion{0.25, material}};
- const Material material_avr = createAveragedMaterial(material, regions);
+ const Material material_avr = CreateAveragedMaterial(material, regions);
EXPECT_EQ(material_avr.getName(), material.getName() + "_avg");
EXPECT_EQ(material_avr.magnetization(), magnetization);
EXPECT_DOUBLE_EQ(material_avr.materialData().real(), 0.5);
@@ -171,7 +171,7 @@ TEST_F(MaterialTest, AveragedMaterialTest)
EXPECT_TRUE(material_avr.typeID() == MATERIAL_TYPES::RefractiveMaterial);
const Material material2 = MaterialBySLD();
- const Material material_avr2 = createAveragedMaterial(material2, regions);
+ const Material material_avr2 = CreateAveragedMaterial(material2, regions);
const complex_t expected_res = std::conj(1.0 - std::sqrt(complex_t(0.5, 0.25)));
EXPECT_DOUBLE_EQ(material_avr2.materialData().real(), expected_res.real());
EXPECT_DOUBLE_EQ(material_avr2.materialData().imag(), expected_res.imag());
@@ -179,12 +179,12 @@ TEST_F(MaterialTest, AveragedMaterialTest)
EXPECT_TRUE(material_avr2.typeID() == MATERIAL_TYPES::RefractiveMaterial);
const Material material3 = MaterialBySLD("Material3", 0.5, 0.5, magnetization);
- EXPECT_THROW(createAveragedMaterial(material3, regions), std::runtime_error);
+ EXPECT_THROW(CreateAveragedMaterial(material3, regions), std::runtime_error);
const Material material4 = HomogeneousMaterial();
const std::vector<HomogeneousRegion> regions2
= {HomogeneousRegion{0.25, material3}, HomogeneousRegion{0.25, material3}};
- const Material material_avr3 = createAveragedMaterial(material4, regions2);
+ const Material material_avr3 = CreateAveragedMaterial(material4, regions2);
EXPECT_DOUBLE_EQ(material_avr3.materialData().real(), 0.25);
EXPECT_DOUBLE_EQ(material_avr3.materialData().imag(), 0.25);
EXPECT_EQ(material_avr3.magnetization(), kvector_t(0.5, 0.0, 0.0));
--
GitLab