Fix problem with slicing core shell particles

Core/Particle/Particle.cpp

@@ -69,6 +69,8 @@ SlicedParticle Particle::createSlicedParticle(ZLimits limits) const
         P_rotation.reset(mP_rotation->clone());
     std::unique_ptr<IFormFactor> P_temp_ff(
                 mP_form_factor->createSlicedFormFactor(limits, *P_rotation, m_position));
+    if (!P_temp_ff)
+        return {};
     std::unique_ptr<FormFactorDecoratorMaterial> P_ff(new FormFactorDecoratorMaterial(*P_temp_ff));
     double volume = P_temp_ff->volume();
     Material transformed_material(

Core/Particle/ParticleCoreShell.cpp

@@ -53,8 +53,6 @@ SlicedParticle ParticleCoreShell::createSlicedParticle(ZLimits limits) const
     P_core->rotate(*P_rotation);
     P_core->translate(m_position);
     auto sliced_core = P_core->createSlicedParticle(limits);
-    if (!sliced_core.mP_slicedff || sliced_core.m_regions.size()!=1)
-        return {};
 
     // shell
     std::unique_ptr<Particle> P_shell(mp_shell->clone());
@@ -64,6 +62,14 @@ SlicedParticle ParticleCoreShell::createSlicedParticle(ZLimits limits) const
     if (!sliced_shell.mP_slicedff)
         return {};
 
+    SlicedParticle result;
+    // if core out of limits, return sliced shell
+    if (!sliced_core.mP_slicedff) {
+        result.mP_slicedff.reset(sliced_shell.mP_slicedff.release());
+        result.m_regions.push_back(sliced_shell.m_regions.back());
+        return result;
+    }
+
     // set core ambient material
     if (sliced_shell.m_regions.size()!=1)
         return {};
@@ -71,7 +77,6 @@ SlicedParticle ParticleCoreShell::createSlicedParticle(ZLimits limits) const
     sliced_core.mP_slicedff->setAmbientMaterial(shell_material);
 
     // construct sliced particle
-    SlicedParticle result;
     sliced_shell.m_regions.back().m_volume -= sliced_core.m_regions.back().m_volume;
     result.mP_slicedff.reset(new FormFactorCoreShell(sliced_core.mP_slicedff.release(),
                                                      sliced_shell.mP_slicedff.release()));

Core/Scattering/IFormFactor.cpp

@@ -25,6 +25,8 @@
 namespace {
 bool ShapeIsContainedInLimits(const IFormFactor& formfactor, ZLimits limits,
                               const IRotation& rot, kvector_t translation);
+bool ShapeOutsideLimits(const IFormFactor& formfactor, ZLimits limits,
+                        const IRotation& rot, kvector_t translation);
 }
 
 IFormFactor::~IFormFactor() {}
@@ -34,6 +36,8 @@ IFormFactor* IFormFactor::createSlicedFormFactor(ZLimits limits, const IRotation
 {
     if (ShapeIsContainedInLimits(*this, limits, rot, translation))
         return CreateTransformedFormFactor(*this, rot, translation);
+    if (ShapeOutsideLimits(*this, limits, rot, translation))
+        return nullptr;
     if (canSliceAnalytically(rot))
         return sliceFormFactor(limits, rot, translation);
     throw std::runtime_error(getName() + "::createSlicedFormFactor error: not supported for "
@@ -96,4 +100,17 @@ bool ShapeIsContainedInLimits(const IFormFactor& formfactor, ZLimits limits,
         return false;
     return true;
 }
+bool ShapeOutsideLimits(const IFormFactor& formfactor, ZLimits limits,
+                        const IRotation& rot, kvector_t translation)
+{
+    double zbottom = formfactor.bottomZ(rot) + translation.z();
+    double ztop = formfactor.topZ(rot) + translation.z();
+    OneSidedLimit lower_limit = limits.lowerLimit();
+    OneSidedLimit upper_limit = limits.upperLimit();
+    if (!upper_limit.m_limitless && zbottom >= upper_limit.m_value)
+        return true;
+    if (!lower_limit.m_limitless && ztop <= lower_limit.m_value)
+        return true;
+    return false;
+}
 }