Skip to content

Instantly share code, notes, and snippets.

@IgorBaratta
Last active February 27, 2024 14:12
Show Gist options
  • Select an option

  • Save IgorBaratta/8fedb06c5000fd2967b3053044461823 to your computer and use it in GitHub Desktop.

Select an option

Save IgorBaratta/8fedb06c5000fd2967b3053044461823 to your computer and use it in GitHub Desktop.
test_assemble_scalar
cmake_minimum_required(VERSION 3.18)
set(PROJECT_NAME test_scalar)
project(${PROJECT_NAME} LANGUAGES C CXX)
# Set C++20 standard
set(CMAKE_CXX_STANDARD 20)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF)
find_package(DOLFINX REQUIRED)
# Add target to compile UFL files
set(SCALAR_TYPE "--scalar_type=float64")
add_custom_command(
OUTPUT form.c
COMMAND ffcx ${CMAKE_CURRENT_SOURCE_DIR}/form.py ${SCALAR_TYPE}
VERBATIM
DEPENDS form.py
COMMENT "Compile form.py using FFCx"
)
set(CMAKE_INCLUDE_CURRENT_DIR ON)
add_executable(${PROJECT_NAME} main.cpp ${CMAKE_CURRENT_BINARY_DIR}/form.c)
target_link_libraries(${PROJECT_NAME} dolfinx)
from basix.ufl import element
from ufl import (
Coefficient,
FunctionSpace,
Mesh,
Measure,
inner,
)
e = element("Lagrange", "triangle", 1)
coord_element = element("Lagrange", "triangle", 1, shape=(2,))
mesh = Mesh(coord_element)
dx = Measure("dx", domain=mesh)
V = FunctionSpace(mesh, e)
f = Coefficient(V)
M = inner(f, f) * dx(1)
#include "form.h"
#include <cmath>
#include <dolfinx.h>
#include <vector>
using namespace dolfinx;
using T = double;
using U = typename dolfinx::scalar_value_type_t<T>;
int main(int argc, char *argv[])
{
dolfinx::init_logging(argc, argv);
MPI_Init(&argc, &argv);
MPI_Comm comm = MPI_COMM_WORLD;
{
auto mesh = std::make_shared<mesh::Mesh<U>>(
mesh::create_rectangle<U>(comm, {{{0.0, 0.0}, {1.0, 1.0}}},
{10, 10}, mesh::CellType::triangle));
auto topology = mesh->topology();
int tdim = topology->dim();
// Create function space
auto V = std::make_shared<fem::FunctionSpace<U>>(
fem::create_functionspace(functionspace_form_form_M, "f", mesh));
std::int32_t num_cells = topology->index_map(tdim)->size_local();
std::vector<std::int32_t> cell_indices(num_cells);
std::iota(cell_indices.begin(), cell_indices.end(), 0);
std::vector<std::int32_t> cell_values(num_cells, 0.0);
for (std::int32_t i = 0; i < num_cells / 2; i++)
cell_values[i] = 1.0;
auto itype = fem::IntegralType::cell;
std::vector<std::pair<int, std::vector<std::int32_t>>> cell_domains = fem::compute_integration_domains(itype, *topology, cell_indices, tdim, cell_values);
std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>> domains;
for (auto &[domain, indices] : cell_domains)
domains.push_back({domain, indices});
// Create a function in the function space
auto f = std::make_shared<fem::Function<U>>(V);
// set values to 1 and compute area
std::fill(f->x()->mutable_array().begin(),
f->x()->mutable_array().end(), 1.0);
// Create form
auto M = std::make_shared<fem::Form<T, U>>(fem::create_form<T, U>(
*form_form_M, {}, {{"f", f}}, {}, {{itype, domains}}, mesh));
// Assemble form
T area = dolfinx::fem::assemble_scalar(*M);
MPI_Allreduce(MPI_IN_PLACE, &area, 1, MPI_DOUBLE, MPI_SUM, comm);
// Print area
if (dolfinx::MPI::rank(comm) == 0)
std::cout << "Area: " << area << std::endl;
}
MPI_Finalize();
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment