Create a SubDomain

[EmilyBourne]

I would like to be able to create some kind of SubDomain. For example this could be achieved by adding a stride to the Domain. This would be useful in several places.

It can be used to iterate over a sliced domain. This can be useful for creating a Lagrange interpolation or calculating a quadrature (e.g. Simpson's quadrature).
E.g. currently a Simpson's quadrature is written as:

    for (auto it = domain.begin() + 1; it < domain.end() - 1; it += 2) {
        ddc::DiscreteElement<IDim> idx = *it;
        coefficients(idx) = 2. / 3. * (distance_at_left(idx) + distance_at_right(idx));
        idx += ddc::DiscreteVector<IDim>(1);
        coefficients(idx) = 1. / 3. * (distance_at_left(idx) + distance_at_right(idx));
    }

with a SubDomain it would be:

    for (ddc::DiscreteElement<IDim> idx : subdomain) {
        ddc::DiscreteElement<IDim> idx2 = idx + ddc::DiscreteVector<IDim>(1);
        coefficients(idx) = 2. / 3. * (distance_at_left(idx) + distance_at_right(idx));
        coefficients(idx2) = 1. / 3. * (distance_at_left(idx2) + distance_at_right(idx2));
    }

It can also be used to create a Chunk on a SubDomain. Eg when working with patches it is common to need values of a function over all the domain, but the derivatives only at certain points (usually the boundary). Currently this is only possible by explicitly choosing the number of points at the compilation, or creating a separate mesh for the derivatives (which then have different DiscreteElements)

[tpadioleau]

I am thinking of defining a ddc::Chunk/ddc::ChunkSpan not only over a ddc::DiscreteDomain but over a cartesian product of 1D ranges of ddc::DiscreteElement. In this case ddc::DiscreteDomain would be the particular case of a cartesian product of 1D continuous ranges. This way we could also introduce a strided range of ddc::DiscreteElement.

Accessing ddc::Chunk/ddc::ChunkSpan data using ddc::DiscreteElement could be slow depending on the type of range. For other reasons I would also like to introduce an access using ddc::DiscreteVector that would always be "fast" compared to the latter.

Do you think it would help you ?

[EmilyBourne]

This sounds like it would cover our use cases as described in the issue. I'm not sure I have fully understood though. The Chunk/ChunkSpan would depend on both a DiscreteDomain and a Range?

I can see a use for 3 kinds of ranges:

1. continuous range (equivalent to the current behaviour)
2. strided range (as described in the issue)
3. arbitrary (ordered) range (general case)

@pauline987654321 has recently come across a use case for the generala case.

For other reasons I would also like to introduce an access using ddc::DiscreteVector that would always be "fast" compared to the latter.

This sounds like a good idea in theory but I wonder if ddc::DiscreteVector is the right tool?
It seems a little confusing to me if you have:

NewRangeType valid_indices(start_element, n_elements, stride_length);
Chunk my_chunk(valid_indices);
ddc::DiscreteVector step(3);
my_chunk(start_element + step) != my_chunk(step) // start_element + step != start_element + step * stride_length

[tpadioleau]

This sounds like it would cover our use cases as described in the issue. I'm not sure I have fully understood though. The Chunk/ChunkSpan would depend on both a DiscreteDomain and a Range?

No I was thinking of only one parameter describing the Support of the Chunk/ChunkSpan, for example something like ChunkSpan<double, StridedDomain<DDimX, DDimY>>.

3. arbitrary (ordered) range (general case)

Yes, it seems more difficult but should be doable. However accessing it using DiscreteElement would be more costly because of the lookup.

This sounds like a good idea in theory but I wonder if ddc::DiscreteVector is the right tool? It seems a little confusing to me if you have:

NewRangeType valid_indices(start_element, n_elements, stride_length);
Chunk my_chunk(valid_indices);
ddc::DiscreteVector step(3);
my_chunk(start_element + step) != my_chunk(step) // start_element + step != start_element + step * stride_length

To me you would have to wonder if the element belongs to the range, which is not the case in your example start_element + step does not belong to valid_indices. Then you cannot call my_chunk(start_element + step) because it does not define any attribute associated to the element start_element + step. The third element in valid_indices would be valid_indices(step)/valid_indices[step] and this element could be used to get the attribute with my_chunk(valid_indices(step)), even though my_chunk(step) would be both easier to read and faster.

[EmilyBourne]

No I was thinking of only one parameter

Sounds good

Yes, it seems more difficult but should be doable. However accessing it using DiscreteElement would be more costly because of the lookup.

Yes that makes sense

To me you would have to wonder if the element belongs to the range, which is not the case in your example start_element + step does not belong to valid_indices

If stride_length is 3 then start_element + step does belong to valid_indices

valid_indices(step)/valid_indices[step]

So you would use operator() for access via indices and operator[] for access via offset? That makes sense but could definitely be a trip hasard for new users

[tpadioleau]

If stride_length is 3 then start_element + step does belong to valid_elements

Oh that is true, then it would be the second element in valid_elements, so something like this should work

ddc::DiscreteVector range_stride(3);
StridedDomain valid_elements(start_element, n_elements, range_stride);
ddc::Chunk<int, StridedDomain> my_chunk(valid_elements);
EXPECT_EQUAL(my_chunk(ddc::DiscreteVector(0)), my_chunk(start_element));
EXPECT_EQUAL(my_chunk(ddc::DiscreteVector(1)), my_chunk(start_element + range_stride));
EXPECT_EQUAL(my_chunk(ddc::DiscreteVector(1)), my_chunk(valid_elements(ddc::DiscreteVector(1))));

Would this make sense ?

valid_elements(step)/valid_elements[step]

So you would use operator() for access via indices and operator[] for access via offset? That makes sense but could definitely be a trip hasard for new users

No I meant one of the two operators would allow to access the n-th element using a DiscreteVector, I just don't know yet which one.

[EmilyBourne]

Then this looks good to me. Presumably an iterator would then use a ddc::DiscreteVector internally?

[tpadioleau]

Could be yes

[EmilyBourne]

Then that all sounds ok to me. What would a for_each loop look like in this case?

ddc::DiscreteVector range_stride(2);
StridedDomain valid_elements(start_element, n_elements, range_stride);
ddc::Chunk<int, StridedDomain> my_chunk(valid_elements);

// Option 1
for_each (valid_elements, [&](ddc::DiscreteElement<IDim> idx) {
     my_chunk(idx) = ddc::coordinate(idx);
}

// Option 2
for_each (valid_elements, [&](ddc::DiscreteVector<IDim> vec) {
     my_chunk(valid_elements(vec)) = ddc::coordinate(start_element + vec);
}

// Option 3
for_each (valid_elements, [&](ChunkIterator iter) {
     ddc::DiscreteElement idx = iter.element();
     my_chunk(iter) = ddc::coordinate(idx);
}

Something else?

[tpadioleau]

I think we can define algorithms for the options 1 and 2.

[tpadioleau]

I can see a use for 3 kinds of ranges:

1. continuous range (equivalent to the current behaviour)
2. strided range (as described in the issue)
3. arbitrary (ordered) range (general case)

@pauline987654321 has recently come across a use case for the generala case.

Back to this comment, in 3. would it require storage or would there be an analytical formula ?

[Paulinemoi]

It would require storage. We won't always have an analytical formula.

[tpadioleau]

@Paulinemoi Hello, as you may have seen we have a first proof of concept in #698. Feel free to have a look at it, especially the tests, and feel free to give us feedback on the feature.

[PaulineVidal]

In the test TEST(StorageDiscreteDomainTest, Constructor)
https://github.com/CExA-project/ddc/pull/698/files#diff-697ad69ab98aff5e69baead89914f4bb45e7af9640e1e81a890abef61e6e185bR86,
I am not sure to get why it is
EXPECT_EQ(dom_xy.distance_from_front(lbound_x + 1, lbound_y + 0), DVectXY(1, 0));
and not
EXPECT_EQ(dom_xy.distance_from_front(lbound_x + 2, lbound_y + 0), DVectXY(1, 0));.

DDomXY const dom_xy(view_x, view_y); is a stride domain with the indices in {0, 2}x{0,2,5}, right?
So, with the definition of view_x, is dom_xy.distance_from_front(lbound_x + 1, lbound_y + 0) working?

[tpadioleau]

Thank you for taking time to look at it.

I am not sure to get why it is EXPECT_EQ(dom_xy.distance_from_front(lbound_x + 1, lbound_y + 0), DVectXY(1, 0)); and not EXPECT_EQ(dom_xy.distance_from_front(lbound_x + 2, lbound_y + 0), DVectXY(1, 0));.

Because it seems there is a typo :) unfortunately I think the tests were passing :(

DDomXY const dom_xy(view_x, view_y); is a stride domain with the indices in {0, 2}x{0,2,5}, right?

Actually the DDomY cannot be represented as a strided domain justifying the new StorageDiscreteDomain. This domain selects the elements lbound_x, lbound_x + 2 and lbound_x + 5. It could be represented as a StridedDiscreteDomain if we replace lbound_x + 5 by lbound_x + 4.

So, with the definition of view_x, is dom_xy.distance_from_front(lbound_x + 1, lbound_y + 0) working?

No, only DiscreteElement<DDimY> contained in dom_xy can be used to read the data from the associated ChunkSpan.