Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

introduce MixedMesh class #303

Open
wants to merge 1 commit into
base: main
Choose a base branch
from

Conversation

ksagiyam
Copy link
Contributor

@ksagiyam ksagiyam commented Aug 21, 2024

First half of #264

Introduce MixedMesh class for multi-mesh problems.

Note that extract_domains() now uses more robust sort_domains() inside, so it might behave slightly differently.

Edit 12-09-2024:

MixedMesh class represents a collection of meshes (e.g., submeshes) that, along with a MixedElement, can represent a mixed function space defined across multiple domains.
The motivation is to allow for treating Arguments and Coefficients on a mixed function space defined across multiple domains just like those on a mixed function space defined on a single mesh.

Specifically, the following becomes possible (see tests for more):

    cell = triangle
    mesh0 = Mesh(FiniteElement("Lagrange", cell, 1, (2,), identity_pullback, H1))
    mesh1 = Mesh(FiniteElement("Lagrange", cell, 1, (2,), identity_pullback, H1))
    domain = MixedMesh([mesh0, mesh1])
    elem0 = FiniteElement("Lagrange", cell, 1, (), identity_pullback, H1)
    elem1 = FiniteElement("Lagrange", cell, 2, (), identity_pullback, H1)
    elem = MixedElement([elem0, elem1])
    V = FunctionSpace(domain, elem)
    v = TestFunction(V)
    v0, v1 = split(v)

For now, one can only perform cell integrations when MixedMeshes are used. This is because, e.g., an interior facet integration on a submesh may either be interior or exterior facet integration on the parent mesh, and we need a logic to apply default restrictions on coefficients defined on the participating meshes. This is the second half of #264.

Also, currently, all component meshes must have the same cell type (and thus the same topological dimension) -- we are to remove this limitation in the future.

Core changes:

  • GradRuleSet.{reference_value, reference_grad} work component-wise (component of the mixed space) if the FunctionSpace is defined on a MixedMesh, so that each component is associated with a component of the MixedMesh, say domain_i (JacobianInverse(domain_i) is then well defined).
  • extract_arguments_and_coefficients is now called extract_terminals_with_domain, and it now also collects GeometricQuantys, so that we can correctly handle, say SpatialCoordinate(mesh1) and SpatialCoordinate(mesh2), in problem solving environments.

ufl/domain.py Outdated Show resolved Hide resolved
@ksagiyam ksagiyam force-pushed the ksagiyam/introduce_mixed_map_0 branch from 0eeb02f to 71d79f4 Compare August 21, 2024 13:16
@ksagiyam
Copy link
Contributor Author

ksagiyam commented Sep 6, 2024

@jpdean @jorgensd This is what we briefly discussed at PDESoft. Could you have a look?

base_coeff_and_args = extract_type(a, (BaseArgument, BaseCoefficient))
arguments = [f for f in base_coeff_and_args if isinstance(f, BaseArgument)]
coefficients = [f for f in base_coeff_and_args if isinstance(f, BaseCoefficient)]
base_coeff_and_args_and_gq = extract_type(a, (BaseArgument, BaseCoefficient, GeometricQuantity))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could we just call this list base_types, as it is quite explicit from the function call what it does, and it makes the latter code more readable.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Renamed it terminals.

Comment on lines +240 to +239
if not isinstance(gq._domain, Mesh):
raise TypeError(f"{gq}._domain must be a Mesh: got {gq._domain}")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What else com gq._domain be? is this relying on the deprecated definition that a domain can be defined just through a ufl.Cell?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No. It checks if gq._domain is Mesh and not MixedMesh.

Comment on lines 34 to 35
domain = MixedMesh(mesh0, mesh1, mesh2)
V = FunctionSpace(domain, elem)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So we now have two ways of creating a FunctionSpace with MixedElements.

  1. A single mesh of class Mesh with a MixedElement, where MixedElement can have arbitrary many elements
  2. A MixedMesh consisting of M meshes, and a MixedElement consisting of M elements.

Do we need additional checks in the FunctionSpace constructor to ensure that this is satisfied (instead of catching this at a later instance, inside for instance GenericDerivativeRuleset.reference_value.

I would also like to highlight that the API for MixedElement and MixedMeshes differs, as one passes a list, while the other unrolls a list, I think we should be consistent, and use lists in both places

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Added checks in FunctionSpace constructor and changed MixedMesh API as suggested. Thanks.

if element.num_sub_elements != len(domain):
raise RuntimeError(f"{element.num_sub_elements} != {len(domain)}")
g = ReferenceGrad(o)
vsh = g.ufl_shape[:-1]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could we give this a more obvious name?
Isn't this isn't this just the ufl_shape of the input operand, which should always be the same as the ufl_shape of o.
Another question (as I don't quite see how this works) is what is the ref_dim?, which would be the topological dimension of the unique "ufl-domain" of o. But what is the topological dimension of a MixedMesh?

Does this mean that a MixedMesh cannot have co-dimension 1 meshes in them, i.e. a Mesh of triangles and a Mesh of intervals?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Turned out vsh could be removed, so I did.

You are right. Currently, MixedMesh can only be composed of meshes of the same cell_type, so topological dimension is well defined. I edited the doc string of MixedMesh to make this point clear.

@garth-wells
Copy link
Member

It would be helpful to have a PR summary, e.g. what it sets out the achieve, what it supports, designs considered, design rationale, implementation summary, and what a MixedMesh is and how a MixedMesh is different from a Mesh.

@ksagiyam ksagiyam force-pushed the ksagiyam/introduce_mixed_map_0 branch from 6dbbd72 to 3512b14 Compare September 12, 2024 12:46
@ksagiyam
Copy link
Contributor Author

Updated PR summary.

Added more checks.

Fixed MixedPullback.physical_value_shape().

@ksagiyam ksagiyam force-pushed the ksagiyam/introduce_mixed_map_0 branch from 3512b14 to 510af2f Compare October 15, 2024 14:03
@ksagiyam
Copy link
Contributor Author

Can I have another round of review on this?

@ksagiyam ksagiyam force-pushed the ksagiyam/introduce_mixed_map_0 branch from 510af2f to 9c1e697 Compare November 15, 2024 22:11
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants