@@ -44,54 +44,71 @@ class GeometricProcessStrategy : public ProcessStrategy<NumericType, D> {
4444 return ProcessResult::INVALID_INPUT ;
4545 }
4646
47- auto maskLevelSet = SmartPointer<viennals::Domain<NumericType, D>>::New (
48- domain->getGrid ());
49- bool foundMaterial = false ;
50-
51- auto const &levelSets = domain->getLevelSets ();
52- for (auto &material : maskMaterials) {
53- for (int j = 0 ; j < levelSets.size (); ++j) {
54- if (materialMap->getMaterialAtIdx (j) == material) {
55- auto lsCopy = SmartPointer<viennals::Domain<NumericType, D>>::New (
56- levelSets[j]);
57-
58- // remove all lower level sets that are not mask materials
59- for (int k = j - 1 ; k >= 0 ; --k) {
60- if (MaterialMap::isMaterial (materialMap->getMaterialAtIdx (k),
61- maskMaterials))
62- continue ;
63-
64- viennals::BooleanOperation<NumericType, D>(
65- lsCopy, levelSets[k],
66- viennals::BooleanOperationEnum::RELATIVE_COMPLEMENT )
67- .apply ();
68- }
69-
70- if (foundMaterial) {
71- // union with mask level set
72- viennals::BooleanOperation<NumericType, D>(
73- maskLevelSet, lsCopy, viennals::BooleanOperationEnum::UNION )
74- .apply ();
75- } else {
76- maskLevelSet = lsCopy;
77- foundMaterial = true ;
47+ if (geometricModel->isSingleMaterial () && maskMaterials.size () == 1 ) {
48+ // operation should be applied to a single material, so we use the
49+ // surface as mask and exclude the material from the surface
50+ mask = viennals::Domain<NumericType, D>::New (domain->getSurface ());
51+ auto materialLS = domain->getMaterialLevelSet (maskMaterials[0 ]);
52+ if (materialLS) {
53+ viennals::BooleanOperation<NumericType, D>(
54+ mask, materialLS,
55+ viennals::BooleanOperationEnum::RELATIVE_COMPLEMENT )
56+ .apply ();
57+ VIENNACORE_LOG_DEBUG (
58+ " Created mask level set from single mask material: " +
59+ materialMap->toString (maskMaterials[0 ]));
60+ }
61+ } else {
62+ auto maskLevelSet =
63+ viennals::Domain<NumericType, D>::New (domain->getGrid ());
64+ bool foundMaterial = false ;
65+ auto const &levelSets = domain->getLevelSets ();
66+ for (auto &material : maskMaterials) {
67+ for (int j = 0 ; j < levelSets.size (); ++j) {
68+ if (materialMap->getMaterialAtIdx (j) == material) {
69+ auto lsCopy =
70+ viennals::Domain<NumericType, D>::New (levelSets[j]);
71+
72+ // remove all lower level sets that are not mask materials
73+ for (int k = j - 1 ; k >= 0 ; --k) {
74+ if (MaterialMap::isMaterial (materialMap->getMaterialAtIdx (k),
75+ maskMaterials))
76+ continue ;
77+
78+ viennals::BooleanOperation<NumericType, D>(
79+ lsCopy, levelSets[k],
80+ viennals::BooleanOperationEnum::RELATIVE_COMPLEMENT )
81+ .apply ();
82+ }
83+
84+ if (foundMaterial) {
85+ // union with mask level set
86+ viennals::BooleanOperation<NumericType, D>(
87+ maskLevelSet, lsCopy,
88+ viennals::BooleanOperationEnum::UNION )
89+ .apply ();
90+ } else {
91+ maskLevelSet = lsCopy;
92+ foundMaterial = true ;
93+ }
7894 }
7995 }
8096 }
81- }
8297
83- if (maskLevelSet->getNumberOfPoints () > 0 ) {
84- mask = maskLevelSet;
85- if (Logger::hasDebug ()) {
86- auto dbgMesh = viennals::Mesh<NumericType>::New ();
87- viennals::ToMesh<NumericType, D>(mask, dbgMesh).apply ();
88- viennals::VTKWriter<NumericType>(dbgMesh, " geometric_mask_debug" )
89- .apply ();
98+ if (maskLevelSet->getNumberOfPoints () > 0 ) {
99+ mask = maskLevelSet;
100+ if (Logger::hasDebug ()) {
101+ auto dbgMesh = viennals::Mesh<NumericType>::New ();
102+ viennals::ToMesh<NumericType, D>(mask, dbgMesh).apply ();
103+ viennals::VTKWriter<NumericType>(dbgMesh, " geometric_mask_debug" )
104+ .apply ();
105+ }
106+ } else {
107+ VIENNACORE_LOG_WARNING (
108+ " None of the specified mask materials were found in the "
109+ " domain, "
110+ " cannot create mask level set from mask materials." );
90111 }
91- } else {
92- VIENNACORE_LOG_WARNING (
93- " None of the specified mask materials were found in the domain, "
94- " cannot create mask level set from mask materials." );
95112 }
96113 }
97114 }
@@ -104,8 +121,7 @@ class GeometricProcessStrategy : public ProcessStrategy<NumericType, D> {
104121 // last level set.
105122 for (int i = context.domain ->getNumberOfLevelSets () - 1 ; i >= 0 ; --i) {
106123 viennals::BooleanOperation<NumericType, D>(
107- context.domain ->getLevelSets ()[i],
108- context.domain ->getLevelSets ().back (),
124+ context.domain ->getLevelSets ()[i], context.domain ->getSurface (),
109125 viennals::BooleanOperationEnum::INTERSECT )
110126 .apply ();
111127 }
0 commit comments