Skip to content

Conversation

@Rohit-Kakodkar
Copy link
Collaborator

Description

Please describe the changes/features in this pull request.

  • Adds a dimension tag to materials
  • Adds element tags within the mesher
  • Adds constructor to element types
  • Adds constructor to assembly properties

Issue Number

If there is an issue created for these changes, link it here

Checklist

Please make sure to check developer documentation on specfem docs.

  • I ran the code through pre-commit to check style
  • THE DOCUMENTATION BUILDS WITHOUT WARNINGS/ERRORS
  • I have added labels to the PR (see right hand side of the PR page)
  • My code passes all the integration tests
  • I have added sufficient unittests to test my changes
  • I have added/updated documentation for the changes I am proposing
  • I have updated CMakeLists to ensure my code builds
  • My code builds across all platforms

- [x] Adds a dimension tag to materials
- [x] Adds element tags within the mesher
- [x] Adds constructor to element types
- [x] Adds constructor to assembly properties
- [x] Adds new constructor for receivers
- [x] MESHFEM3D mesher now writes number of NGLL points
- [x] Added element grid to mesh struct
- [x] Added test to check element grid
@Rohit-Kakodkar Rohit-Kakodkar added the enhancement New feature or request label Oct 30, 2025
@lsawade lsawade linked an issue Oct 31, 2025 that may be closed by this pull request
4 tasks
*/
template <specfem::element::medium_tag MediumTag>
struct material<MediumTag, specfem::element::property_tag::anisotropic>
template <specfem::dimension::type DimensionTag,
Copy link
Collaborator

Choose a reason for hiding this comment

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

2 and 3D stiffness are different!
https://en.wikipedia.org/wiki/Voigt_notation

In 2D we omit storing all the ones that are related to out-of-plane strain.

From specfem3d

c11 = c11store(INDEX_IJK,ispec)
        c12 = c12store(INDEX_IJK,ispec)
        c13 = c13store(INDEX_IJK,ispec)
        c14 = c14store(INDEX_IJK,ispec)
        c15 = c15store(INDEX_IJK,ispec)
        c16 = c16store(INDEX_IJK,ispec)
        c22 = c22store(INDEX_IJK,ispec)
        c23 = c23store(INDEX_IJK,ispec)
        c24 = c24store(INDEX_IJK,ispec)
        c25 = c25store(INDEX_IJK,ispec)
        c26 = c26store(INDEX_IJK,ispec)
        c33 = c33store(INDEX_IJK,ispec)
        c34 = c34store(INDEX_IJK,ispec)
        c35 = c35store(INDEX_IJK,ispec)
        c36 = c36store(INDEX_IJK,ispec)
        c44 = c44store(INDEX_IJK,ispec)
        c45 = c45store(INDEX_IJK,ispec)
        c46 = c46store(INDEX_IJK,ispec)
        c55 = c55store(INDEX_IJK,ispec)
        c56 = c56store(INDEX_IJK,ispec)
        c66 = c66store(INDEX_IJK,ispec)

        sigma_xx = c11 * duxdxl(INDEX_IJK) + c16 * duxdyl_plus_duydxl + c12 * duydyl(INDEX_IJK) + &
                   c15 * duzdxl_plus_duxdzl + c14 * duzdyl_plus_duydzl + c13 * duzdzl(INDEX_IJK)
        sigma_yy = c12 * duxdxl(INDEX_IJK) + c26 * duxdyl_plus_duydxl + c22 * duydyl(INDEX_IJK) + &
                   c25 * duzdxl_plus_duxdzl + c24 * duzdyl_plus_duydzl + c23 * duzdzl(INDEX_IJK)
        sigma_zz = c13 * duxdxl(INDEX_IJK) + c36 * duxdyl_plus_duydxl + c23 * duydyl(INDEX_IJK) + &
                   c35 * duzdxl_plus_duxdzl + c34 * duzdyl_plus_duydzl + c33 * duzdzl(INDEX_IJK)
        sigma_xy = c16 * duxdxl(INDEX_IJK) + c66 * duxdyl_plus_duydxl + c26 * duydyl(INDEX_IJK) + &
                   c56 * duzdxl_plus_duxdzl + c46 * duzdyl_plus_duydzl + c36 * duzdzl(INDEX_IJK)
        sigma_xz = c15 * duxdxl(INDEX_IJK) + c56 * duxdyl_plus_duydxl + c25 * duydyl(INDEX_IJK) + &
                   c55 * duzdxl_plus_duxdzl + c45 * duzdyl_plus_duydzl + c35 * duzdzl(INDEX_IJK)
        sigma_yz = c14 * duxdxl(INDEX_IJK) + c46 * duxdyl_plus_duydxl + c24 * duydyl(INDEX_IJK) + &
                   c45 * duzdxl_plus_duxdzl + c44 * duzdyl_plus_duydzl + c34 * duzdzl(INDEX_IJK)

struct material {
int n_materials; ///< Number of elements
std::vector<specfem::medium::material<type, property> >
std::vector<specfem::medium::material<dimension, type, property> >
Copy link
Collaborator

Choose a reason for hiding this comment

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

Suggested change
std::vector<specfem::medium::material<dimension, type, property> >
std::vector<specfem::medium::material<dimension_tag, type, property> >

Needs updating above and below

specfem::mesh::materials<specfem::dimension::type::dim2>::material<type, property>::material(
const int n_materials,
const std::vector<specfem::medium::material<type, property> >
const std::vector<specfem::medium::material<dimension, type, property> >
Copy link
Collaborator

Choose a reason for hiding this comment

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

Suggested change
const std::vector<specfem::medium::material<dimension, type, property> >
const std::vector<specfem::medium::material<dimension_tag, type, property> >

* @return std::variant Material properties
*/
specfem::medium::material<MediumTag, PropertyTag>
specfem::medium::material<dimension, MediumTag, PropertyTag>
Copy link
Collaborator

Choose a reason for hiding this comment

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

Suggested change
specfem::medium::material<dimension, MediumTag, PropertyTag>
specfem::medium::material<dimension_tag, MediumTag, PropertyTag>

* Each dimension defines the number of GLL points used for numerical
* integration and interpolation in that coordinate direction.
*/
struct GLLGrid {
Copy link
Collaborator

Choose a reason for hiding this comment

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

element_grid<dim3>?

Comment on lines 117 to 199
FOR_EACH_IN_PRODUCT(
(DIMENSION_TAG(DIM3), MEDIUM_TAG(ELASTIC), PROPERTY_TAG(ISOTROPIC),
BOUNDARY_TAG(NONE)),
CAPTURE(elements, h_elements) {
int count = 0;
int index = 0;
for (int ispec = 0; ispec < nspec; ispec++) {
if (medium_tags(ispec) == _medium_tag_) {
count++;
}
}
_elements_ =
IndexViewType("specfem::assembly::element_types::elements", count);
_h_elements_ = Kokkos::create_mirror_view(_elements_);
for (int ispec = 0; ispec < nspec; ispec++) {
if (medium_tags(ispec) == _medium_tag_) {
_h_elements_(index) = ispec;
index++;
}
}
Kokkos::deep_copy(_elements_, _h_elements_);
})

FOR_EACH_IN_PRODUCT(
(DIMENSION_TAG(DIM3), MEDIUM_TAG(ELASTIC), PROPERTY_TAG(ISOTROPIC),
BOUNDARY_TAG(NONE)),
CAPTURE(elements, h_elements) {
int count = 0;
int index = 0;

for (int ispec = 0; ispec < nspec; ispec++) {
if (medium_tags(ispec) == _medium_tag_ &&
property_tags(ispec) == _property_tag_) {
count++;
}
}

_elements_ =
IndexViewType("specfem::assembly::element_types::elements", count);
_h_elements_ = Kokkos::create_mirror_view(_elements_);

for (int ispec = 0; ispec < nspec; ispec++) {
if (medium_tags(ispec) == _medium_tag_ &&
property_tags(ispec) == _property_tag_) {
_h_elements_(index) = ispec;
index++;
}
}

Kokkos::deep_copy(_elements_, _h_elements_);
})

FOR_EACH_IN_PRODUCT(
(DIMENSION_TAG(DIM3), MEDIUM_TAG(ELASTIC), PROPERTY_TAG(ISOTROPIC),
BOUNDARY_TAG(NONE)),
CAPTURE(elements, h_elements) {
int count = 0;
int index = 0;

for (int ispec = 0; ispec < nspec; ispec++) {
if (medium_tags(ispec) == _medium_tag_ &&
property_tags(ispec) == _property_tag_ &&
boundary_tags(ispec) == _boundary_tag_) {
count++;
}
}

_elements_ =
IndexViewType("specfem::assembly::element_types::elements", count);
_h_elements_ = Kokkos::create_mirror_view(_elements_);

for (int ispec = 0; ispec < nspec; ispec++) {
if (medium_tags(ispec) == _medium_tag_ &&
property_tags(ispec) == _property_tag_ &&
boundary_tags(ispec) == _boundary_tag_) {
_h_elements_(index) = ispec;
index++;
}
}

Kokkos::deep_copy(_elements_, _h_elements_);
})
}
Copy link
Collaborator

Choose a reason for hiding this comment

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

This seems redundant? Something is not right, you are setting elements_dim3_elastic_isotropic_none 3 times with different counts, overwriting the first two. Am I missing somethign here?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Good catch. I had copied the function from above which has the same bug.

Copy link
Collaborator

@icui icui left a comment

Choose a reason for hiding this comment

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

Minor formatting, otherwise look good to me.

*/
bool operator==(
const material<MediumTag,
const material<dimension_tag, MediumTag,
Copy link
Collaborator

Choose a reason for hiding this comment

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

Suggested change
const material<dimension_tag, MediumTag,
const material<dimension_tag, medium_tag,

bool operator==(
const material<MediumTag, specfem::element::property_tag::isotropic>
&other) const {
const material<dimension_tag, MediumTag,
Copy link
Collaborator

Choose a reason for hiding this comment

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

Suggested change
const material<dimension_tag, MediumTag,
const material<dimension_tag, medium_tag,

* @tparam PropertyTag Property tag for the material
*/
template <specfem::element::medium_tag MediumTag,
template <specfem::dimension::type dimension_tag,
Copy link
Collaborator

Choose a reason for hiding this comment

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

Suggested change
template <specfem::dimension::type dimension_tag,
template <specfem::dimension::type DimensionTag,

@Rohit-Kakodkar Rohit-Kakodkar merged commit 51be740 into devel Nov 7, 2025
8 checks passed
@Rohit-Kakodkar Rohit-Kakodkar deleted the issue-1281 branch November 7, 2025 18:46
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

enhancement New feature or request

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Compute material properties using MESHFEM3D mesh

4 participants