Hi all,
Just wanted to follow up on this in case it could help someone else looking
to do something similar.
Based on Winni's very helpful response, I was able to successfully create a
mesh of the unit square
$[0,1]^2$ that contains the interface described by $\Gamma(x) = (x, 0.5 +
0.1 \sin(2\pi x))$.
Please see the attached--it's not a standalone .cpp file (i.e. not a
'minimal working example'), but, it
should be able to be copied directly into Step-49 and work right away.
Thanks again--
--Sean
On Tuesday, September 10, 2024 at 5:52:09 PM UTC-4 Sean Carney wrote:
> Dear Winni,
>
> Thank you very much for pointing me to the helpful example from the code
> from your recent paper.
>
> I have learned a great deal from breaking down this function and
> understanding what each part is
> doing (thank you for adding the descriptive comments within).
>
> I think I now understand reasonably well how to create a mesh with some
> kind of internal interface $\Gamma$
> *so long as* the interface is described by one of the built-in functions
> within the GridGenerator namespace
> (for example, hyper_ball, cylinder_shell, etc).
>
> In this case one can use the built-in Manifold objects, such as
> SphericalManifold (as you use), PolarManifold, etc.
>
> However, I'm still not sure how to create an internal interface with my
> own specified function. E.g. something
> as simple as creating the domain [-1,1]x[0,1] with the interface $\Gamma$
> parameterized by points along the
> curve $(0.2 \sin(2\pi y), y)$.
>
> I guess the FunctionManifold
> <https://www.dealii.org/current/doxygen/deal.II/classFunctionManifold.html>
> class is my best bet, but I'm still trying to understand how this works.
>
> If you have any further suggestions, they would be appreciated.
> But thank you again for the helpful comment.
>
> Best regards,
> -Sean
>
> On Friday, September 6, 2024 at 2:59:04 AM UTC-4 winnifried wrote:
>
>> Dear Sean,
>> I think you can find an example for what you want to achieve in the code
>> of a recent paper https://epubs.siam.org/doi/10.1137/23M158262X
>>
>> The code is available at
>>
>> https://zenodo.org/records/10459309
>>
>> where you can check the file
>> InverseHomogenization/geometry.h
>>
>> Best regards
>> Winni
>>
>> On 05.09.24 22:13, Sean Carney wrote:
>> > Hi all,
>> >
>> > I want to create a 2D domain in deal.II with an internal, curve
>> interface
>> > $\Gamma$. A simple example of what I want to create is in the attached
>> > screenshot (credit to Bartels, Bonito & Tscherner
>> > <https://doi.org/10.1093/imanum/drad004>).
>> >
>> > In short, does anyone have any suggestions for how to do this?
>> > The critical thing is that upon mesh refinement, deal.II knows where my
>> > interface is and adds points accordingly (of course, this issue is
>> already
>> > discussed in Step-1). Here, however, the 1D manifold is internal, so
>> maybe
>> > this changes things?
>> >
>> > What I tried so far is the following:
>> >
>> > I have an analytic expression for my curve (i.e. it can be written in
>> the
>> > form $(\Gamma(y), y)$), so I created two different triangulation
>> objects
>> > using the
>> > GridGenerator::subdivided_hyper_rectangle
>> > and
>> > GridTools::transform
>> > functions. Then, I tried to merge them together. However, I am
>> receiving an
>> > error:
>> >
>> > "An error occurred in line <1235> of file
>> > <dealii/source/grid/grid_tools.cc> in function
>> > void dealii::GridTools::{anonymous}::AdjacentCells<2>::push_back(const
>> > dealii::GridTools::{anonymous}::AdjacentCell&)
>> > The violated condition was:
>> > adjacent_cells[1].cell_index == numbers::invalid_unsigned_int"
>> >
>> > Of course, I know that the vertices along the common interface, in this
>> > case $\Gamma$, need to be aligned. I believe that this is indeed the
>> case.
>> >
>> > My code is simple enough:
>> >
>> >
>> > void grid_10()
>> > {
>> > const int N = 14;
>> >
>> > // create 1st triang
>> > Triangulation<2> tri_1;
>> >
>> > std::vector<unsigned int> repetitions(2);
>> > repetitions[0] = N;
>> > repetitions[1] = N;
>> > GridGenerator::subdivided_hyper_rectangle(tri_1,
>> > repetitions,
>> > Point<2>(0.0, 0.0),
>> > Point<2>(2./3., 1.));
>> >
>> > GridTools::transform(
>> > [](const Point<2> &in) {
>> > if (in[0] < 0.5)
>> > {
>> > return in;
>> > }
>> > else
>> > {
>> > return Point<2>(2./3. - 1./6.*std::sin(numbers::PI * in[1]), in[1]);
>> > }
>> > },
>> > tri_1);
>> >
>> > // create 2nd triang
>> > Triangulation<2> tri_2;
>> >
>> > repetitions[0] = N;
>> > repetitions[1] = N;
>> > GridGenerator::subdivided_hyper_rectangle(tri_2,
>> > repetitions,
>> > Point<2>(2./3., 0),
>> > Point<2>(1., 1.));
>> >
>> > GridTools::transform(
>> > [](const Point<2> &in) {
>> > if (in[0] > 0.7)
>> > {
>> > return in;
>> > }
>> > else
>> > {
>> > return Point<2>(2./3. - 1./6.*std::sin(numbers::PI * in[1]), in[1]);
>> > }
>> > },
>> > tri_2);
>> >
>> > // merge them
>> > Triangulation<2> triangulation;
>> > GridGenerator::merge_triangulations(tri_1,tri_2,triangulation);
>> > print_mesh_info(triangulation, "grid.vtu");
>> >
>> > }
>> >
>> > Is there an obvious error that I am missing?
>> >
>> > Even more critically, will this even accomplish what I want? I see in
>> the
>> > "merge_triangulations" documentation that " In particular, manifolds
>> > attached to the two input triangulations will be lost in the result
>> > triangulation."
>> >
>> > Any help is greatly appreciated. Thank you!
>> > Best regards,
>> > -Sean
>> >
>>
>> --
>> --------------------------------------------------------
>> Prof. Dr. Winnifried Wollner
>> Universität Hamburg
>> MIN Faculty
>> Department of Mathematics
>> Bundesstr. 55; 20146 Hamburg; Germany
>> --------------------------------------------------------
>> Room : 115
>> Phone: +49 (0)40 42838-4079 <+49%2040%20428384079>
>> Web : https://www.math.uni-hamburg.de/personen/wollner/
>> eMail: [email protected]
>> --------------------------------------------------------
>>
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/3a29d313-f7ec-4280-897b-0b7fb5d07d57n%40googlegroups.com.
void grid_12()
{
/* create cells from vertices */
const std::vector<Point<2>> vertices{
{0.0, 0.0}, // 0
{0.0, 0.5}, // 1
{0.0, 1.0}, // 2
{0.5, 0.0}, // 3
{0.5, 0.5}, // 4
{0.5, 1.0}, // 5
{1.0, 0.0}, // 6
{1.0, 0.5}, // 7
{1.0, 1.0}, // 8
};
std::vector<CellData<2>> cells(4);
{
cells[0].vertices = {0, 3, 1, 4};
cells[1].vertices = {1, 4, 2, 5};
cells[2].vertices = {3, 6, 4, 7};
cells[3].vertices = {4, 7, 5, 8};
}
Triangulation<2> triangulation;
triangulation.create_triangulation(vertices, cells, SubCellData());
/* Colorize boundaries: */
for (auto face : triangulation.active_face_iterators()) {
const auto center = face->center();
constexpr double eps = 1.0e-6;
if (center[0] < eps) {
face->set_boundary_id(0);
} else if (center[0] > 1.0 - eps) {
face->set_boundary_id(2);
} else if (center[1] < eps) {
face->set_boundary_id(1);
} else if (center[1] > 1.0 - eps) {
face->set_boundary_id(3);
}
}
/*
* Attach
* - a flat manifold on the outer boundaries,
* - a FunctionManifold to faces on the interface,
* - a transfinite interpolation manifold to the rest.
*/
triangulation.set_all_manifold_ids(2);
triangulation.set_all_manifold_ids_on_boundary(numbers::flat_manifold_id);
for (auto face : triangulation.active_face_iterators())
{
double eps = 1.e-6;
const auto center = face->center();
double distance = center[1] - 0.5;
if (std::abs(distance) < eps)
{
face->set_manifold_id(1);
}
}
FunctionManifold<2,2,1> my_func_manifold("x; 0.5 + 0.1*sin(2.*3.14159*x)", "x");
triangulation.set_manifold(1, my_func_manifold);
TransfiniteInterpolationManifold<2> transfinite;
transfinite.initialize(triangulation);
triangulation.set_manifold(2, transfinite);
triangulation.refine_global(3);
print_mesh_info(triangulation, "grid_12.vtu"); // function from Step-49
}
