diff --git a/dev/Requirements.txt b/dev/Requirements.txt new file mode 100644 index 000000000..98f82d1e0 --- /dev/null +++ b/dev/Requirements.txt @@ -0,0 +1,141 @@ +aiofiles==22.1.0 +aiosqlite==0.19.0 +anyio==3.6.2 +appnope==0.1.3 +argon2-cffi==21.3.0 +argon2-cffi-bindings==21.2.0 +arrow==1.2.3 +asttokens==2.2.1 +attrs==23.1.0 +Babel==2.12.1 +backcall==0.2.0 +beautifulsoup4==4.12.2 +bleach==6.0.0 +blinker==1.4 +certifi==2022.12.7 +cffi==1.15.1 +charset-normalizer==3.1.0 +cmake==3.24.1.1 +cmake-generators==1.0.9 +colorama==0.4.4 +comm==0.1.3 +commonmark==0.9.1 +cycler==0.11.0 +debugpy==1.6.7 +decorator==5.1.1 +defusedxml==0.7.1 +distro==1.7.0 +docutils==0.18.1 +executing==1.2.0 +fastjsonschema==2.16.3 +feedgenerator==2.0.0 +fonttools==4.29.1 +fqdn==1.5.1 +future-fstrings==1.2.0 +gcovr==5.2 +gitdb==4.0.9 +GitPython==3.1.27 +h5py==3.8.0 +idna==3.4 +importlib-metadata==6.6.0 +importlib-resources==5.12.0 +ipykernel==6.22.0 +ipython==8.12.1 +ipython-genutils==0.2.0 +isoduration==20.11.0 +jedi==0.18.2 +Jinja2==3.0.3 +joblib==1.3.1 +json5==0.9.11 +jsonpointer==2.3 +jsonschema==4.17.3 +jupyter-events==0.6.3 +jupyter-ydoc==0.2.4 +jupyter_client==8.2.0 +jupyter_core==5.3.0 +jupyter_server==2.5.0 +jupyter_server_fileid==0.9.0 +jupyter_server_terminals==0.4.4 +jupyter_server_ydoc==0.8.0 +jupyterlab==3.6.3 +jupyterlab-pygments==0.2.2 +jupyterlab_server==2.22.1 +kiwisolver==1.3.2 +lxml==4.9.2 +MarkupSafe==2.0.1 +matplotlib==3.5.1 +matplotlib-inline==0.1.6 +mistune==2.0.5 +mpi4py==3.1.4 +nbclassic==0.5.6 +nbclient==0.7.4 +nbconvert==7.3.1 +nbformat==5.8.0 +nest-asyncio==1.5.6 +networkx==3.1 +ninja==1.10.2.4 +nose==1.3.7 +notebook==6.5.4 +notebook_shim==0.2.3 +numpy==1.22.2 +packaging==21.3 +pandas==1.4.1 +pandocfilters==1.5.0 +parso==0.8.3 +pelican==4.7.2 +pexpect==4.8.0 +pickleshare==0.7.5 +Pillow==9.0.1 +pkgutil_resolve_name==1.3.10 +platformdirs==3.5.0 +prometheus-client==0.16.0 +prompt-toolkit==3.0.38 +psutil==5.9.5 +ptyprocess==0.7.0 +pure-eval==0.2.2 +pybind11==2.10.0 +pycparser==2.21 +Pygments==2.11.2 +pyparsing==3.0.7 +pyphen==0.13.2 +pyrsistent==0.19.3 +python-dateutil==2.8.2 +python-json-logger==2.0.7 +pytz==2021.3 +PyYAML==6.0 +pyzmq==25.0.2 +requests==2.29.0 +rfc3339-validator==0.1.4 +rfc3986-validator==0.1.1 +rich==11.2.0 +rodin @ file:///Users/carlos/Projects/rodin +scikit-build==0.15.0 +scikit-learn==1.3.0 +scipy==1.9.0 +seaborn==0.12.2 +Send2Trash==1.8.2 +six==1.16.0 +sklearn==0.0.post7 +smmap==5.0.0 +sniffio==1.3.0 +soupsieve==2.4.1 +stack-data==0.6.2 +svn==1.0.1 +termcolor==2.3.0 +terminado==0.17.1 +threadpoolctl==3.2.0 +tinycss2==1.2.1 +tomli==2.0.1 +tornado==6.3.1 +traitlets==5.9.0 +typing_extensions==4.5.0 +Unidecode==1.3.2 +uri-template==1.2.0 +urllib3==1.26.15 +wcwidth==0.2.6 +webcolors==1.13 +webencodings==0.5.1 +websocket-client==1.5.1 +y-py==0.5.9 +ypy-websocket==0.8.2 +zipp==3.15.0 diff --git a/dev/py/include_dependencies.py b/dev/py/include_dependencies.py new file mode 100644 index 000000000..9e8894296 --- /dev/null +++ b/dev/py/include_dependencies.py @@ -0,0 +1,241 @@ +import os +import re +import numpy as np +import matplotlib.pyplot as plt +import argparse +import networkx as nx +from pathlib import Path +from shutil import which +from termcolor import colored +from networkx.drawing.nx_agraph import write_dot + +PT_TO_INCH = (1.0 / 72) + +MIN_NODE_SIZE = 2.0 +MAX_NODE_SIZE = 7.0 + +MIN_FONT_SIZE = 6 +MAX_FONT_SIZE = 42 + + +COMMON_CPP_EXTENSIONS = ['.h', '.hpp', '.hxx'] + +COLORSCHEME = 'spectral9' +NUMBER_OF_COLORS = 9 + +RODIN_TOP_LEVEL_INCLUDES = { + 'Rodin.h', + 'Rodin/Alert.h', + 'Rodin/Variational.h', + 'Rodin/Geometry.h', + 'RodinExternal.h', + 'RodinExternal/MMG.h', +} + +RODIN_TOP_LEVEL_INCLUDES_COLOR_MAP = { + 'Alert.h': 'red', + 'Variational.h': 'green', + 'Geometry.h': 'blue' +} + + +class GraphvizNotFoundError(Exception): + pass + + +def merge_matching_paths(path1, path2): + parts1 = Path(path1).parts + parts2 = Path(path2).parts + if parts2[:len(parts1)] == parts1: + merged_parts = parts1 + parts2[len(parts1):] + merged_path = Path(*merged_parts) + return str(merged_path) + else: + joined_path = Path(path1) / path2 + return str(joined_path) + + +def extract_includes(file_path, context_path): + includes = [] + with open(file_path, 'r') as file: + for line in file: + # system_include_match = re.match(r'#include\s*[<](.*)[>]', line) + user_include_match = re.match(r'#include\s*["](.*)["]', line) + if user_include_match: + include = user_include_match.group(1) + if include == 'ForwardDecls.h': + continue + include_path = os.path.dirname(include) + if include_path: + if not include_path.startswith('Rodin'): + include = merge_matching_paths(context_path, include) + else: + include = os.path.join(context_path, include) + includes.append(include) + return includes + + +def generate_dependency_graph(root_path, extensions): + G = nx.DiGraph() + for root, _, files in os.walk(root_path): + for file in files: + if any(file.endswith(ext) for ext in extensions): + source_path = os.path.join(root, file) + relative_path = os.path.relpath(source_path, root_path) + context_path = os.path.dirname(relative_path) + includes = extract_includes(source_path, context_path) + G.add_node(relative_path) + for include in includes: + G.add_edge(relative_path, include) + return G + +def get_predecessors_map(graph): + sinks = set() + node_to_set = {} + for node in graph.nodes(): + node_to_set[node] = set() + if graph.out_degree(node) == 0: + sinks.add(node) + for sink in sinks: + visited = set() + q = { sink } + while len(q) > 0: + node = q.pop() + visited.add(node) + for n in graph.predecessors(node): + node_to_set[n].add(node) + node_to_set[n] |= node_to_set[node] + if n not in visited: + q.add(n) + return node_to_set + +def get_accumulated_map(graph): + predecessors_map = get_predecessors_map(graph) + node_to_accumulated = {} + for node, predecessors in predecessors_map.items(): + node_to_accumulated[node] = len(predecessors) + return node_to_accumulated + +def get_flow_map(graph): + node_to_flow = {} + for node in graph.nodes(): + out_degree = graph.out_degree(node) + in_degree = graph.in_degree(node) + node_to_flow[node] = in_degree - out_degree + return node_to_flow + +def get_laplacian_map(graph): + laplacian_matrix = nx.directed_laplacian_matrix(graph) + node_to_laplacian = {} + for node_idx, node in enumerate(graph.nodes()): + laplacian_value = laplacian_matrix[node_idx, node_idx] + node_to_laplacian[node] = laplacian_value + return node_to_laplacian + +def get_pagerank_map(graph): + pagerank_scores = nx.pagerank(graph, alpha=0.85) + return pagerank_scores + +def get_color_map(m): + max_v = max(m.values()) + min_v = min(m.values()) + color_map = {} + for node, v in m.items(): + std = float(v - min_v) / float(max_v - min_v) + color_map[node] = int(NUMBER_OF_COLORS - std * (NUMBER_OF_COLORS - 1)) + return color_map + +def get_nodesize_map(m): + max_v = max(m.values()) + min_v = min(m.values()) + size_map = {} + for node, v in m.items(): + std = float(v - min_v) / float(max_v - min_v) + size_map[node] = int((MAX_NODE_SIZE - MIN_NODE_SIZE) * std + MIN_NODE_SIZE) + return size_map + +def get_fontsize_map(m): + max_v = max(m.values()) + min_v = min(m.values()) + fontsize_map = {} + for node, v in m.items(): + std = float(v - min_v) / float(max_v - min_v) + fontsize_map[node] = int((MAX_FONT_SIZE - MIN_FONT_SIZE) * std + MIN_FONT_SIZE) + return fontsize_map + + +def check_graphviz_command(command): + dot_command = which(command) + if not dot_command: + raise GraphvizNotFoundError( + "Graphviz %s command not found. Please install Graphviz to generate graphical output." % command) + +if __name__ == '__main__': + parser = argparse.ArgumentParser( + description='Generates #include dependency graph.') + parser.add_argument('root_path', help='Root path of the codebase.') + parser.add_argument('--layout_engine', + default='dot', + help='Layout engine.') + parser.add_argument('--dot', + help='Output DOT file path.') + parser.add_argument('--png', + help='Output PNG file path.') + parser.add_argument('--svg', + help='Output SVG file path.') + parser.add_argument('--extensions', nargs='+', + default=COMMON_CPP_EXTENSIONS, + help='List of file extensions to consider.') + args = parser.parse_args() + + root_path = args.root_path + folder_name = os.path.basename(os.path.normpath(root_path)) + + dot_filename = args.dot if args.dot else f'{folder_name}_include_dependency_graph.dot' + + try: + check_graphviz_command(args.layout_engine) + except GraphvizNotFoundError as e: + print(colored(e, 'red')) + exit(1) + + root_path = args.root_path + dependency_graph = generate_dependency_graph(root_path, args.extensions) + weight_map = get_accumulated_map(dependency_graph) + color_map = get_color_map(weight_map) + fontsize_map = get_fontsize_map(weight_map) + nodesize_map = get_nodesize_map(weight_map) + + nx.set_node_attributes(dependency_graph, 'filled', name='style') + nx.set_node_attributes(dependency_graph, 'circle', name='shape') + nx.set_node_attributes(dependency_graph, 'true', name='fixedsize') + nx.set_node_attributes(dependency_graph, COLORSCHEME, name='colorscheme') + nx.set_node_attributes(dependency_graph, color_map, name='fillcolor') + nx.set_node_attributes(dependency_graph, weight_map, name='weight') + nx.set_node_attributes(dependency_graph, nodesize_map, name='width') + nx.set_node_attributes(dependency_graph, nodesize_map, name='height') + nx.set_node_attributes(dependency_graph, fontsize_map, name='fontsize') + + write_dot(dependency_graph, dot_filename) + print(colored(f'DOT output generated: {dot_filename}', 'green')) + + layout_engine = args.layout_engine + print(colored(f'Using layout engine {layout_engine}', 'blue')) + command_line = f'{layout_engine} ' + command_line += '-Goverlap=prism -Gbeautify=true -Gsmoothing=true' + if args.png: + png_filename = args.png + command_line += " -Tpng -o %s " % png_filename + + if args.svg: + svg_filename = args.svg + command_line += " -Tsvg -o %s " % svg_filename + + command_line += dot_filename + os.system(command_line) + + if args.png: + print(colored(f'PNG output generated: {png_filename}', 'green')) + if args.svg: + print(colored(f'SVG output generated: {svg_filename}', 'green')) + diff --git a/examples/Misc/RS2023/Domain.cpp b/examples/Misc/RS2023/Domain.cpp index 0ab1ab4ee..bc1d913ba 100644 --- a/examples/Misc/RS2023/Domain.cpp +++ b/examples/Misc/RS2023/Domain.cpp @@ -18,7 +18,7 @@ static constexpr Attribute interior = 1; static constexpr Attribute exterior = 2; static constexpr Attribute ball = 3; -static constexpr Scalar hmax = 0.005; +static constexpr Scalar hmax = 0.01; static constexpr Scalar hmin = 0.01 * hmax; static const Math::Vector x0{{0.5, 0.5}}; // Center of domain diff --git a/examples/Misc/RS2023/Helmholtz/Simulation.cpp b/examples/Misc/RS2023/Helmholtz/Simulation.cpp index 8508ddadd..414eac282 100644 --- a/examples/Misc/RS2023/Helmholtz/Simulation.cpp +++ b/examples/Misc/RS2023/Helmholtz/Simulation.cpp @@ -16,7 +16,7 @@ using namespace Rodin::Math; using namespace Rodin::Geometry; using namespace Rodin::Variational; -static constexpr Scalar hmax = 0.005; // Maximal size of a triangle's edge +static constexpr Scalar hmax = 0.01; // Maximal size of a triangle's edge static constexpr Scalar hmin = 0.1 * hmax; static constexpr Attribute dQ = 2; // Attribute of box boundary @@ -64,20 +64,20 @@ void run(int i, const std::vector& grid); int main(int, char**) { // Define evaluation grid - Math::Vector m_r = Math::Vector::LinSpaced(51, 0, 50); - Math::Vector epsilon_r = Math::Vector::LinSpaced(100, hmax, 0.2); + Math::Vector m_r = Math::Vector::LinSpaced(300, 0, 100); + Math::Vector epsilon_r = Math::Vector::LinSpaced(50, hmax, 0.2); // Math::Vector waveNumber_r = Math::Vector::LinSpaced(1.0 / hmax, 1, 1.0 / hmax); // Math::Vector conductivity_r{{ 1e-12, 0.5, 1.0, 2.0, 1e12 }}; - Math::Vector waveNumber_r = Math::Vector::LinSpaced(81, 0, 80); - // Math::Vector waveNumber_r{{ 0, 1, 2, 3, 4, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 60, 80, 90 }}; + // Math::Vector waveNumber_r = Math::Vector::LinSpaced(150, 0, 74); + Math::Vector waveNumber_r{{ 1, 2, 3, 4, 5, 10, 11, 12, 13, 14, 15, 20 }}; Math::Vector conductivity_r{{ 2.0 }}; Math::Vector angle_r{{ M_PI / 4 }}; std::vector grid; grid.reserve(epsilon_r.size() * waveNumber_r.size() * conductivity_r.size()); - for (const Scalar m : m_r) - for (const Scalar epsilon : epsilon_r) - for (const Scalar waveNumber : waveNumber_r) + for (const Scalar waveNumber : waveNumber_r) + for (const Scalar m : m_r) + for (const Scalar epsilon : epsilon_r) for (const Scalar conductivity : conductivity_r) for (const Scalar angle : angle_r) grid.push_back({ m, epsilon, waveNumber, conductivity, angle }); @@ -167,7 +167,11 @@ void run(int id, const std::vector& grid) P1 vh(mesh); // Define oscillatory screen - ScalarFunction h = 2; + ScalarFunction h = + [&](const Point& p) + { + return 2 + sin(2 * pi * data.m * p.x()) * sin(2 * pi * data.m * p.y()); + }; // Define conductivity ScalarFunction gamma = 1; diff --git a/examples/Misc/RS2023/L2_Grid_Plot.py b/examples/Misc/RS2023/L2_Grid_Plot.py index c3e52c3db..a061c710f 100644 --- a/examples/Misc/RS2023/L2_Grid_Plot.py +++ b/examples/Misc/RS2023/L2_Grid_Plot.py @@ -9,11 +9,11 @@ plt.style.use("bmh") conductivities = [2.0] -waveNumbers = [ 0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 60, 80, 90 ] +waveNumbers = [1, 2, 3, 4, 5, 10, 11, 12, 13, 14, 15, 20] angles = [ math.pi / 4 ] -x_column = 'waveNumber' +x_column = 'm' y_column = 'epsilon' latex_label = { @@ -38,46 +38,53 @@ df = pd.read_csv(args.filename) print('Read dataset !') - for angle in angles: - for conductivity in conductivities: - plt.figure() - sdf = df[ - (df['conductivity'] > conductivity - 0.001) & - (df['conductivity'] < conductivity + 0.001) & - (df['angle'] > angle - 0.001) & - (df['angle'] < angle + 0.001) & - (df['waveNumber'] < 10) - ] + thresh = math.exp(df['error'].median()) - thresh = sdf['error'].median() - sdf.loc[sdf['error'] > thresh, 'error'] = thresh + df.loc[df['error'] > thresh, 'error'] = thresh - plt.hexbin( - sdf.loc[:, x_column], - sdf.loc[:, y_column], - gridsize=25, - reduce_C_function=np.mean, - cmap='inferno', - C=sdf.loc[:, 'error'], edgecolors='face', linewidths=0.5) + i = 0 + for waveNumber in waveNumbers: + print(waveNumber) + for angle in angles: + for conductivity in conductivities: + plt.figure() + sdf = df[ + (df['conductivity'] > conductivity - 0.001) & + (df['conductivity'] < conductivity + 0.001) & + (df['angle'] > angle - 0.001) & + (df['angle'] < angle + 0.001) & + (df['waveNumber'] > waveNumber - 0.01) & + (df['waveNumber'] < waveNumber + 0.01) + ] - # plt.contourf( - # sdf.loc[:, 'm'], - # sdf.loc[:, 'epsilon'], - # sdf.loc[:, 'error']) - plt.colorbar() + plt.hexbin( + sdf.loc[:, x_column], + sdf.loc[:, y_column], + gridsize=50, + reduce_C_function=np.mean, + cmap='inferno', + C=sdf.loc[:, 'error'], edgecolors='face', linewidths=0.5) - plt.xlabel(latex_label[x_column]) - plt.ylabel(latex_label[y_column]) - plt.title( - '$|| u_\epsilon - u_0 ||_{L^2 (Q)} \ - \ \mathrm{with} \ \gamma |_{B(x, e)} = %.2E, \ - \ \\theta = %d^\circ $' % (conductivity, (180.0 / math.pi) * angle), - pad=20) - out=("%s_VS_%s_Gamma=%.2E_Angle=%.2E") % ( - file_label[x_column], file_label[y_column], conductivity, (180.0 / math.pi) * angle) - plt.savefig(out + '.svg') - plt.savefig(out + '.png') - # plt.show() + # plt.contourf( + # sdf.loc[:, 'm'], + # sdf.loc[:, 'epsilon'], + # sdf.loc[:, 'error']) + plt.colorbar() + + plt.xlabel(latex_label[x_column]) + plt.ylabel(latex_label[y_column]) + plt.title( + '$|| u_\epsilon - u_0 ||_{L^2 (Q / B(x, 0.25))} \ + \ \mathrm{with} \ \gamma |_{B(x, e)} = %.2E, \ + \ \\theta = %d^\circ, \ k=%.2E$' % ( + conductivity, (180.0 / math.pi) * angle, waveNumber), pad=20) + out=("%s_VS_%s_Gamma=%.2E_Angle=%.2E_Wavenumber=%.2E") % ( + file_label[x_column], file_label[y_column], conductivity, + (180.0 / math.pi) * angle, waveNumber) + plt.savefig(str(i) + '_' + out + '.svg') + plt.savefig(str(i) + '_' + out + '.png') + # plt.show() + i += 1 diff --git a/resources/py/plot.py b/resources/py/plot.py deleted file mode 100644 index 7011858e8..000000000 --- a/resources/py/plot.py +++ /dev/null @@ -1,45 +0,0 @@ -""" -Description ------------ - -Parses and plots the contents of a one column file. This is useful for -generating uniform plots of the files generated by some of the shape -optimization examples. - -Usage ------ -python plot.py - -""" - -import numpy as np -from matplotlib import pyplot as plt -plt.style.use('grayscale') -import numpy as np -import pandas as pd -import argparse - -if __name__ == '__main__': - description = '''Parses and plots the contents of a one column file. This is - useful for generating uniform plots of the files generated by some of the shape - optimization examples.''' - parser = argparse.ArgumentParser(description=description) - parser.add_argument('filename', type=str, - help='Name of file to plot.') - args = parser.parse_args() - - data = pd.read_csv(args.filename) - data = pd.DataFrame(data) - - x, y = data.iloc[1:, 0], data.iloc[1:, 1] - m, b = np.polyfit(np.log10(x), np.log10(y), 1) - print('m: %f, b: %f' % (m, b)) - plt.plot(x, y) - plt.grid(alpha=0.1, aa=True) - plt.xscale('log') - plt.yscale('log') - plt.xlabel(data.columns[0], fontsize=18) - plt.ylabel(data.columns[1], fontsize=18) - plt.savefig(args.filename + '.svg', format='svg') - plt.show() - diff --git a/src/Rodin/Alert/MemberFunctionException.h b/src/Rodin/Alert/MemberFunctionException.h index 51367d083..314acdb32 100644 --- a/src/Rodin/Alert/MemberFunctionException.h +++ b/src/Rodin/Alert/MemberFunctionException.h @@ -8,6 +8,7 @@ #define RODIN_ALERT_MEMBERFUNCTIONEXCEPTION_H #include +#include #include "Exception.h" #include "Identifier.h" diff --git a/src/Rodin/Geometry/ForwardDecls.h b/src/Rodin/Geometry/ForwardDecls.h index 3eec56ed1..ef7599ea2 100644 --- a/src/Rodin/Geometry/ForwardDecls.h +++ b/src/Rodin/Geometry/ForwardDecls.h @@ -56,6 +56,8 @@ namespace Rodin::Geometry template class Mesh; + class SubMeshBase; + /** * @brief Represents a subset of a Mesh. * @@ -80,6 +82,9 @@ namespace Rodin::Geometry template class SubMesh; + template <> + class SubMesh; + template class SubMeshBuilder; } diff --git a/src/Rodin/Geometry/Mesh.cpp b/src/Rodin/Geometry/Mesh.cpp index 1b482eaa1..5418b8d81 100644 --- a/src/Rodin/Geometry/Mesh.cpp +++ b/src/Rodin/Geometry/Mesh.cpp @@ -4,8 +4,7 @@ * (See accompanying file LICENSE or copy at * https://www.boost.org/LICENSE_1_0.txt) */ -#include -#include "Rodin/Alert.h" +#include "Rodin/Alert/MemberFunctionException.h" #include "Rodin/Variational/P1.h" #include "Rodin/Variational/GridFunction.h" @@ -30,46 +29,173 @@ namespace Rodin::Geometry } // ---- Mesh ---------------------------------------------- + Mesh& + Mesh::load(const boost::filesystem::path& filename, IO::FileFormat fmt) + { + std::ifstream input(filename.c_str()); + if (!input) + { + Alert::MemberFunctionException(*this, __func__) + << "Failed to open " << filename << " for reading." + << Alert::Raise; + } + switch (fmt) + { + case IO::FileFormat::MFEM: + { + IO::MeshLoader loader(*this); + loader.load(input); + break; + } + case IO::FileFormat::MEDIT: + { + IO::MeshLoader loader(*this); + loader.load(input); + break; + } + default: + { + Alert::MemberFunctionException(*this, __func__) + << "Loading from \"" << fmt << "\" format unsupported." + << Alert::Raise; + break; + } + } + return *this; + } - Mesh Mesh::UniformGrid(Polytope::Type g, size_t h, size_t w) + void Mesh::save( + const boost::filesystem::path& filename, + IO::FileFormat fmt, size_t precision) const { - Builder build; - const size_t dim = Polytope::getGeometryDimension(g); - build.initialize(dim).nodes(h * w); - switch (g) + std::ofstream ofs(filename.c_str()); + if (!ofs) { - case Polytope::Type::Point: + Alert::MemberFunctionException(*this, __func__) + << "Failed to open " << filename << " for writing." + << Alert::Raise; + } + ofs.precision(precision); + switch (fmt) + { + case IO::FileFormat::MFEM: { - return build.nodes(1).vertex({0}).finalize(); + IO::MeshPrinter printer(*this); + printer.print(ofs); + break; } - case Polytope::Type::Triangle: + case IO::FileFormat::MEDIT: { - assert(h * w >= 4); - for (size_t i = 0; i < h; i++) - { - for (size_t j = 0; j < w; j++) - build.vertex({ static_cast(j), static_cast(i) }); - } + IO::MeshPrinter printer(*this); + printer.print(ofs); + break; + } + default: + { + Alert::MemberFunctionException(*this, __func__) + << "Saving to \"" << fmt << "\" format unsupported." + << Alert::Raise; + } + } + ofs.close(); + } - build.reserve(dim, 2 * (h - 1) * (w - 1)); - for (size_t i = 0; i < h - 1; i++) + SubMesh Mesh::keep(Attribute attr) const + { + return keep(FlatSet{attr}); + } + + SubMesh Mesh::keep(const FlatSet& attrs) const + { + const size_t D = getDimension(); + SubMesh::Builder build; + build.initialize(*this); + for (Index i = 0; i < getElementCount(); i++) + { + if (attrs.count(getAttribute(D, i))) + { + build.include(D, i); + for (size_t d = 1; d <= D - 1; d++) { - for (size_t j = 0; j < w - 1; j++) - { - build.polytope(g, { i * w + j, i * w + j + 1, (i + 1) * w + j }) - .polytope(g, { i * w + j + 1, (i + 1) * w + j + 1, (i + 1) * w + j }); - } + const auto& inc = getConnectivity().getIncidence(D, d); + if (inc.size() > 0) + build.include(d, inc.at(i)); } - return build.finalize(); } - default: + } + return build.finalize(); + } + + SubMesh Mesh::skin() const + { + const size_t D = getDimension(); + SubMesh::Builder build; + build.initialize(*this); + for (auto it = getBoundary(); !it.end(); ++it) + { + const Index i = it->getIndex(); + build.include(D - 1, i); + for (size_t d = 1; d <= D - 2; d++) { - assert(false); - return build.nodes(0).finalize(); + const auto& inc = getConnectivity().getIncidence(D - 1, d); + if (inc.size() > 0) + build.include(d, inc.at(i)); } - }; + } + return build.finalize(); } + SubMesh Mesh::trim(Attribute attr) const + { + return trim(FlatSet{attr}); + } + + SubMesh Mesh::trim(const FlatSet& attrs) const + { + const size_t D = getDimension(); + SubMesh::Builder build; + build.initialize(*this); + for (Index i = 0; i < getElementCount(); i++) + { + if (!attrs.count(getAttribute(D, i))) + { + build.include(D, i); + for (size_t d = 1; d <= D - 1; d++) + { + const auto& inc = getConnectivity().getIncidence(D, d); + if (inc.size() > 0) + build.include(d, inc.at(i)); + } + } + } + return build.finalize(); + } + + Mesh& Mesh::displace(const + Variational::GridFunction>>& u) + { + assert(u.getFiniteElementSpace().getVectorDimension() == getSpaceDimension()); + m_vertices += u.getData(); + flush(); + return *this; + } + + Mesh& Mesh::scale(Scalar c) + { + m_vertices *= c; + flush(); + return *this; + } + +#ifdef RODIN_USE_MPI + Mesh + Mesh::parallelize(boost::mpi::communicator comm) + { + return Mesh(comm, *this); + } +#endif + Eigen::Map Mesh::getVertexCoordinates(Index idx) const { const auto size = static_cast(getSpaceDimension()); @@ -137,49 +263,6 @@ namespace Rodin::Geometry } } - Mesh& Mesh::scale(Scalar c) - { - m_vertices *= c; - flush(); - return *this; - } - - void Mesh::save( - const boost::filesystem::path& filename, - IO::FileFormat fmt, size_t precision) const - { - std::ofstream ofs(filename.c_str()); - if (!ofs) - { - Alert::Exception() - << "Failed to open " << filename << " for writing." - << Alert::Raise; - } - ofs.precision(precision); - switch (fmt) - { - case IO::FileFormat::MFEM: - { - IO::MeshPrinter printer(*this); - printer.print(ofs); - break; - } - case IO::FileFormat::MEDIT: - { - IO::MeshPrinter printer(*this); - printer.print(ofs); - break; - } - default: - { - Alert::Exception() - << "Saving to \"" << fmt << "\" format unsupported." - << Alert::Raise; - } - } - ofs.close(); - } - Scalar MeshBase::getVolume() { Scalar totalVolume = 0; @@ -256,15 +339,6 @@ namespace Rodin::Geometry // return res; // } - // ---- Mesh ------------------------------------------------------ -#ifdef RODIN_USE_MPI - Mesh - Mesh::parallelize(boost::mpi::communicator comm) - { - return Mesh(comm, *this); - } -#endif - size_t Mesh::getCount(size_t dimension) const { return m_connectivity.getCount(dimension); @@ -377,120 +451,69 @@ namespace Rodin::Geometry return *this; } - Mesh& - Mesh::load(const boost::filesystem::path& filename, IO::FileFormat fmt) + SubMeshBase& Mesh::asSubMesh() { - std::ifstream input(filename.c_str()); - if (!input) + assert(isSubMesh()); + if (!isSubMesh()) { - Alert::Exception() - << "Failed to open " << filename << " for reading." + Alert::MemberFunctionException(*this, __func__) + << "This instance of Mesh is not a SubMesh (isSubMesh() == false). " + << "Downcasting to SubMesh is ill-defined." << Alert::Raise; } - switch (fmt) - { - case IO::FileFormat::MFEM: - { - IO::MeshLoader loader(*this); - loader.load(input); - break; - } - case IO::FileFormat::MEDIT: - { - IO::MeshLoader loader(*this); - loader.load(input); - break; - } - default: - { - Alert::Exception() << "Loading from \"" << fmt << "\" format unsupported." - << Alert::Raise; - break; - } - } - - return *this; + return static_cast&>(*this); } - SubMesh Mesh::keep(Attribute attr) const - { - return keep(FlatSet{attr}); - } - - SubMesh Mesh::keep(const FlatSet& attrs) const + const SubMeshBase& Mesh::asSubMesh() const { - const size_t D = getDimension(); - SubMesh::Builder build; - build.initialize(*this); - for (Index i = 0; i < getElementCount(); i++) + assert(isSubMesh()); + if (!isSubMesh()) { - if (attrs.count(getAttribute(D, i))) - { - build.include(D, i); - for (size_t d = 1; d <= D - 1; d++) - { - const auto& inc = getConnectivity().getIncidence(D, d); - if (inc.size() > 0) - build.include(d, inc.at(i)); - } - } + Alert::MemberFunctionException(*this, __func__) + << "This instance of Mesh is not a SubMesh (isSubMesh() == false). " + << "Downcasting to SubMesh is ill-defined." + << Alert::Raise; } - return build.finalize(); + return static_cast&>(*this); } - SubMesh Mesh::skin() const + Mesh Mesh::UniformGrid(Polytope::Type g, size_t h, size_t w) { - const size_t D = getDimension(); - SubMesh::Builder build; - build.initialize(*this); - for (auto it = getBoundary(); !it.end(); ++it) + Builder build; + const size_t dim = Polytope::getGeometryDimension(g); + build.initialize(dim).nodes(h * w); + switch (g) { - const Index i = it->getIndex(); - build.include(D - 1, i); - for (size_t d = 1; d <= D - 2; d++) + case Polytope::Type::Point: { - const auto& inc = getConnectivity().getIncidence(D - 1, d); - if (inc.size() > 0) - build.include(d, inc.at(i)); + return build.nodes(1).vertex({0}).finalize(); } - } - return build.finalize(); - } - - SubMesh Mesh::trim(Attribute attr) const - { - return trim(FlatSet{attr}); - } - - SubMesh Mesh::trim(const FlatSet& attrs) const - { - const size_t D = getDimension(); - SubMesh::Builder build; - build.initialize(*this); - for (Index i = 0; i < getElementCount(); i++) - { - if (!attrs.count(getAttribute(D, i))) + case Polytope::Type::Triangle: { - build.include(D, i); - for (size_t d = 1; d <= D - 1; d++) + assert(h * w >= 4); + for (size_t i = 0; i < h; i++) { - const auto& inc = getConnectivity().getIncidence(D, d); - if (inc.size() > 0) - build.include(d, inc.at(i)); + for (size_t j = 0; j < w; j++) + build.vertex({ static_cast(j), static_cast(i) }); } - } - } - return build.finalize(); - } - Mesh& Mesh::displace(const - Variational::GridFunction>>& u) - { - assert(u.getFiniteElementSpace().getVectorDimension() == getSpaceDimension()); - m_vertices += u.getData(); - flush(); - return *this; + build.reserve(dim, 2 * (h - 1) * (w - 1)); + for (size_t i = 0; i < h - 1; i++) + { + for (size_t j = 0; j < w - 1; j++) + { + build.polytope(g, { i * w + j, i * w + j + 1, (i + 1) * w + j }) + .polytope(g, { i * w + j + 1, (i + 1) * w + j + 1, (i + 1) * w + j }); + } + } + return build.finalize(); + } + default: + { + assert(false); + return build.nodes(0).finalize(); + } + }; } } diff --git a/src/Rodin/Geometry/Mesh.h b/src/Rodin/Geometry/Mesh.h index 24248b290..1aaf15504 100644 --- a/src/Rodin/Geometry/Mesh.h +++ b/src/Rodin/Geometry/Mesh.h @@ -39,8 +39,32 @@ namespace Rodin::Geometry public: virtual ~MeshBase() = default; + inline + constexpr + bool operator==(const MeshBase& other) const + { + return this == &other; + } + + inline + constexpr + bool operator!=(const MeshBase& other) const + { + return this != &other; + } + virtual MeshBase& scale(Scalar c) = 0; + virtual MeshBase& load( + const boost::filesystem::path& filename, + IO::FileFormat fmt = IO::FileFormat::MFEM) = 0; + + virtual void save( + const boost::filesystem::path& filename, + IO::FileFormat fmt = IO::FileFormat::MFEM, size_t precison = 16) const = 0; + + virtual void flush() = 0; + /** * @brief Indicates whether the mesh is a surface or not. * @@ -52,6 +76,44 @@ namespace Rodin::Geometry */ bool isSurface() const; + /** + * @brief Indicates whether the mesh is a submesh or not. + * @returns True if mesh is a submesh, false otherwise. + * + * A Mesh which is also a SubMesh may be casted into down to access + * the SubMesh functionality. For example: + * @code{.cpp} + * if (mesh.isSubMesh()) + * { + * // Cast is well defined + * auto& submesh = static_cast(mesh); + * } + * @endcode + * + */ + virtual bool isSubMesh() const = 0; + + virtual bool isInterface(Index faceIdx) const = 0; + + /** + * @brief Determines whether a face of the mesh is on the boundary. + */ + virtual bool isBoundary(Index faceIdx) const = 0; + + /** + * @brief Gets the dimension of the elements. + * @returns Dimension of the elements. + * @see getSpaceDimension() const + */ + virtual size_t getDimension() const = 0; + + /** + * @brief Gets the dimension of the ambient space + * @returns Dimension of the space which the mesh is embedded in + * @see getDimension() const + */ + virtual size_t getSpaceDimension() const = 0; + /** * @brief Gets the number of vertices in the mesh. */ @@ -130,62 +192,6 @@ namespace Rodin::Geometry */ virtual const FlatSet& getAttributes(size_t d) const = 0; - bool operator==(const MeshBase& other) const - { - return this == &other; - } - - bool operator!=(const MeshBase& other) const - { - return this != &other; - } - - /** - * @brief Gets the dimension of the elements. - * @returns Dimension of the elements. - * @see getSpaceDimension() const - */ - virtual size_t getDimension() const = 0; - - /** - * @brief Gets the dimension of the ambient space - * @returns Dimension of the space which the mesh is embedded in - * @see getDimension() const - */ - virtual size_t getSpaceDimension() const = 0; - - virtual MeshBase& load( - const boost::filesystem::path& filename, - IO::FileFormat fmt = IO::FileFormat::MFEM) = 0; - - virtual void save( - const boost::filesystem::path& filename, - IO::FileFormat fmt = IO::FileFormat::MFEM, size_t precison = 16) const = 0; - - /** - * @brief Indicates whether the mesh is a submesh or not. - * @returns True if mesh is a submesh, false otherwise. - * - * A Mesh which is also a SubMesh may be casted into down to access - * the SubMesh functionality. For example: - * @code{.cpp} - * if (mesh.isSubMesh()) - * { - * // Cast is well defined - * auto& submesh = static_cast(mesh); - * } - * @endcode - * - */ - virtual bool isSubMesh() const = 0; - - virtual bool isInterface(Index faceIdx) const = 0; - - /** - * @brief Determines whether a face of the mesh is on the boundary. - */ - virtual bool isBoundary(Index faceIdx) const = 0; - virtual FaceIterator getBoundary() const = 0; virtual FaceIterator getInterface() const = 0; @@ -220,7 +226,9 @@ namespace Rodin::Geometry virtual Eigen::Map getVertexCoordinates(Index idx) const = 0; - virtual void flush() = 0; + virtual SubMeshBase& asSubMesh() = 0; + + virtual const SubMeshBase& asSubMesh() const = 0; }; /// Index containing the indices of boundary elements. @@ -552,6 +560,10 @@ namespace Rodin::Geometry */ virtual SubMesh keep(const FlatSet& attrs) const; + SubMeshBase& asSubMesh() override; + + const SubMeshBase& asSubMesh() const override; + /** * @brief Loads a mesh from file in the given format. * @param[in] filename Name of file to read diff --git a/src/Rodin/Geometry/SubMesh.cpp b/src/Rodin/Geometry/SubMesh.cpp index 7826d47df..ff2f1170f 100644 --- a/src/Rodin/Geometry/SubMesh.cpp +++ b/src/Rodin/Geometry/SubMesh.cpp @@ -27,7 +27,7 @@ namespace Rodin::Geometry return m_parent.get(); } - Point SubMesh::inclusion(const Point& p) + Point SubMesh::inclusion(const Point& p) const { const auto& parent = getParent(); const auto& polytope = p.getPolytope(); @@ -39,7 +39,7 @@ namespace Rodin::Geometry std::cref(p.getReferenceCoordinates()), p.getPhysicalCoordinates()); } - Point SubMesh::restriction(const Point& p) + Point SubMesh::restriction(const Point& p) const { const auto& child = *this; const auto& polytope = p.getPolytope(); diff --git a/src/Rodin/Geometry/SubMesh.h b/src/Rodin/Geometry/SubMesh.h index 86850969c..2254bb3e8 100644 --- a/src/Rodin/Geometry/SubMesh.h +++ b/src/Rodin/Geometry/SubMesh.h @@ -17,18 +17,33 @@ namespace Rodin::Geometry { - /** - * @defgroup DotSpecializations Dot Template Specializations - * @brief Template specializations of the Dot class. - * @see Dot + * @defgroup SubMeshSpecializations SubMesh Template Specializations + * @brief Template specializations of the SubMesh class. + * @see SubMesh */ - /** - * @ingroup DotSpecializations - */ + class SubMeshBase + { + public: + virtual Point inclusion(const Point& p) const = 0; + + virtual Point restriction(const Point& p) const = 0; + + /** + * @returns Reference to the parent Mesh object + */ + virtual const MeshBase& getParent() const = 0; + + /** + * @brief Gets the map of polytope indices from the SubMesh to the parent + * Mesh. + */ + virtual const boost::bimap& getPolytopeMap(size_t d) const = 0; + }; /** + * @ingroup SubMeshSpecializations * @brief A SubMesh object represents a subregion of a Mesh object. * * A SubMesh object contains a reference to the parent Mesh object. It also @@ -41,13 +56,23 @@ namespace Rodin::Geometry * if (mesh.isSubMesh()) * { * // Cast is well defined - * auto& submesh = static_cast(mesh); + * auto& submesh = static_cast&>(mesh); + * } + * @endcode + * + * Alternatively, you can use the asSubMesh() to access the SubMeshBase + * interface: + * @code{.cpp} + * if (mesh.isSubMesh()) + * { + * // Cast is well defined + * auto& submesh = mesh.asSubMesh(); * } * @endcode * */ template <> - class SubMesh : public Mesh + class SubMesh final : public SubMeshBase, public Mesh { public: using Parent = Mesh; @@ -95,9 +120,9 @@ namespace Rodin::Geometry return *this; } - Point inclusion(const Point& p); + Point inclusion(const Point& p) const override; - Point restriction(const Point& p); + Point restriction(const Point& p) const override; bool isSubMesh() const override { @@ -107,13 +132,14 @@ namespace Rodin::Geometry /** * @returns Reference to the parent Mesh object */ - const Mesh& getParent() const; + const Mesh& getParent() const override; /** * @brief Gets the map of polytope indices from the SubMesh to the parent * Mesh. */ - const boost::bimap& getPolytopeMap(size_t d) const + inline + const boost::bimap& getPolytopeMap(size_t d) const override { return m_s2ps.at(d); } diff --git a/src/Rodin/Variational/FiniteElementSpace.h b/src/Rodin/Variational/FiniteElementSpace.h index 0ec5cc99e..441b13ca8 100644 --- a/src/Rodin/Variational/FiniteElementSpace.h +++ b/src/Rodin/Variational/FiniteElementSpace.h @@ -32,6 +32,20 @@ namespace Rodin::Variational FiniteElementSpaceBase& operator=(FiniteElementSpaceBase&&) = default; + inline + constexpr + bool operator==(const FiniteElementSpaceBase& other) const + { + return this == &other; + } + + inline + constexpr + bool operator!=(const FiniteElementSpaceBase& other) const + { + return this != &other; + } + /** * @brief Gets the total number of degrees of freedom. * @returns Size of the finite element space @@ -110,13 +124,6 @@ namespace Rodin::Variational return static_cast(*this).getMapping(p, v); } }; - - inline - constexpr - bool operator==(const FiniteElementSpaceBase& lhs, const FiniteElementSpaceBase& rhs) - { - return &lhs == &rhs; - } } #endif diff --git a/src/Rodin/Variational/GridFunction.h b/src/Rodin/Variational/GridFunction.h index eaebdbf11..a46facab1 100644 --- a/src/Rodin/Variational/GridFunction.h +++ b/src/Rodin/Variational/GridFunction.h @@ -32,6 +32,17 @@ #include "MatrixFunction.h" #include "FiniteElementSpace.h" +namespace Rodin::FormLanguage +{ + template + struct Traits> + { + using FES = FESType; + using Element = typename Traits::Element; + using RangeType = typename Traits::RangeType; + }; +} + namespace Rodin::Variational { /** @@ -56,18 +67,20 @@ namespace Rodin::Variational * assert(data.cols() == gf.getFiniteElementSpace().getSize()); * ``` */ - template - class GridFunctionBase : public LazyEvaluator> + template ::FES> + class GridFunctionBase : public LazyEvaluator> { public: + using FES = FESType; + /// Type of finite element - using Element = typename FES::Element; + using Element = typename FormLanguage::Traits::Element; /// Parent class using Parent = LazyEvaluator>; /// Range type of value - using RangeType = typename FES::RangeType; + using RangeType = typename FormLanguage::Traits::RangeType; static_assert(std::is_same_v || std::is_same_v); diff --git a/src/Rodin/Variational/P1/GridFunction.h b/src/Rodin/Variational/P1/GridFunction.h index 5ab51c7d8..8695ab2c5 100644 --- a/src/Rodin/Variational/P1/GridFunction.h +++ b/src/Rodin/Variational/P1/GridFunction.h @@ -8,30 +8,49 @@ #define RODIN_VARIATIONAL_P1_GRIDFUNCTION_H #include "Rodin/Variational/GridFunction.h" +#include "Rodin/Geometry/SubMesh.h" + #include "P1.h" +namespace Rodin::FormLanguage +{ + /** + * @ingroup TraitsSpecializations + */ + template + struct Traits>> + { + using RangeType = Range; + using FES = Variational::P1; + using Context = typename FormLanguage::Traits::Context; + using Mesh = typename FormLanguage::Traits::Mesh; + using Element = typename FormLanguage::Traits::Element; + }; +} + namespace Rodin::Variational { /** * @ingroup GridFunctionSpecializations * @brief P1 GridFunction */ - template - class GridFunction> final : public GridFunctionBase>, P1> + template + class GridFunction> final + : public GridFunctionBase>> { public: /// Type of finite element space to which the GridFunction belongs to - using FES = P1; + using FES = P1; /// Type of finite element using Element = typename FES::Element; - /// Parent class - using Parent = GridFunctionBase>, P1>; - /// Range type of value using RangeType = typename FES::RangeType; + /// Parent class + using Parent = GridFunctionBase>; + using Parent::getValue; using Parent::getValueByReference; using Parent::operator=; @@ -78,58 +97,9 @@ namespace Rodin::Variational GridFunction& operator=(const GridFunction&) = delete; - inline - auto getValue(const Geometry::Point& p) const - { - const auto& fes = this->getFiniteElementSpace(); - const auto& polytope = p.getPolytope(); - const size_t d = polytope.getDimension(); - const Index i = polytope.getIndex(); - const auto& fe = fes.getFiniteElement(d, i); - const auto& r = p.getCoordinates(Geometry::Point::Coordinates::Reference); - if constexpr (std::is_same_v) - { - Scalar res = 0; - for (Index local = 0; local < fe.getCount(); local++) - { - const auto& basis = fe.getBasis(local); - res += getValue({d, i}, local) * basis(r); - } - return res; - } - else if constexpr (std::is_same_v) - { - Math::Vector res; - getValueByReference(res, p); - return res; - } - else - { - assert(false); - return void(); - } - } + auto getValue(const Geometry::Point& p) const; - inline - void getValueByReference(Math::Vector& res, const Geometry::Point& p) const - { - static_assert(std::is_same_v); - const auto& fes = this->getFiniteElementSpace(); - const auto& polytope = p.getPolytope(); - const size_t d = polytope.getDimension(); - const Index i = polytope.getIndex(); - const auto& fe = fes.getFiniteElement(d, i); - const auto& r = p.getCoordinates(Geometry::Point::Coordinates::Reference); - const size_t vdim = fes.getVectorDimension(); - const size_t dofs = fe.getCount(); - res.resize(vdim); - res.setZero(); - for (Index local = 0; local < dofs; local++) - { - const auto& basis = fe.getBasis(local); - res += getValue({d, i}, local).coeff(local % vdim) * basis(r); - } - } + void getValueByReference(Math::Vector& res, const Geometry::Point& p) const; inline GridFunction& setWeights() @@ -196,6 +166,8 @@ namespace Rodin::Variational template GridFunction(const P1&) -> GridFunction>; - } + +#include "GridFunction.hpp" + #endif diff --git a/src/Rodin/Variational/P1/GridFunction.hpp b/src/Rodin/Variational/P1/GridFunction.hpp new file mode 100644 index 000000000..a511ff77b --- /dev/null +++ b/src/Rodin/Variational/P1/GridFunction.hpp @@ -0,0 +1,66 @@ +#ifndef RODIN_VARIATIONAL_P1_GRIDFUNCTION_HPP +#define RODIN_VARIATIONAL_P1_GRIDFUNCTION_HPP + +#include "GridFunction.h" + +namespace Rodin::Variational +{ + template + auto GridFunction> + ::getValue(const Geometry::Point& p) const + { + const auto& fes = this->getFiniteElementSpace(); + const auto& fesMesh = fes.getMesh(); + const auto& polytope = p.getPolytope(); + assert(fesMesh == polytope.getMesh()); + const size_t d = polytope.getDimension(); + const Index i = polytope.getIndex(); + const auto& fe = fes.getFiniteElement(d, i); + const auto& r = p.getCoordinates(Geometry::Point::Coordinates::Reference); + if constexpr (std::is_same_v) + { + Scalar res = 0; + for (Index local = 0; local < fe.getCount(); local++) + { + const auto& basis = fe.getBasis(local); + res += getValue({d, i}, local) * basis(r); + } + return res; + } + else if constexpr (std::is_same_v) + { + Math::Vector res; + getValueByReference(res, p); + return res; + } + else + { + assert(false); + return void(); + } + } + + template + void GridFunction> + ::getValueByReference(Math::Vector& res, const Geometry::Point& p) const + { + static_assert(std::is_same_v); + const auto& fes = this->getFiniteElementSpace(); + const auto& polytope = p.getPolytope(); + const size_t d = polytope.getDimension(); + const Index i = polytope.getIndex(); + const auto& fe = fes.getFiniteElement(d, i); + const auto& r = p.getCoordinates(Geometry::Point::Coordinates::Reference); + const size_t vdim = fes.getVectorDimension(); + const size_t dofs = fe.getCount(); + res.resize(vdim); + res.setZero(); + for (Index local = 0; local < dofs; local++) + { + const auto& basis = fe.getBasis(local); + res += getValue({d, i}, local).coeff(local % vdim) * basis(r); + } + } +} + +#endif diff --git a/src/Rodin/Variational/P1/P1.h b/src/Rodin/Variational/P1/P1.h index 01133afdf..80ab7f7bd 100644 --- a/src/Rodin/Variational/P1/P1.h +++ b/src/Rodin/Variational/P1/P1.h @@ -20,6 +20,27 @@ #include "ForwardDecls.h" #include "P1Element.h" +namespace Rodin::FormLanguage +{ + template + struct Traits> + { + using RangeType = Scalar; + using Context = ContextType; + using Mesh = MeshType; + using Element = Variational::ScalarP1Element; + }; + + template + struct Traits> + { + using RangeType = Math::Vector; + using Context = ContextType; + using Mesh = MeshType; + using Element = Variational::VectorP1Element; + }; +} + namespace Rodin::Variational { /** diff --git a/src/Rodin/Variational/P1/P1Element.h b/src/Rodin/Variational/P1/P1Element.h index 030d775af..4819692d7 100644 --- a/src/Rodin/Variational/P1/P1Element.h +++ b/src/Rodin/Variational/P1/P1Element.h @@ -30,6 +30,18 @@ #include "ForwardDecls.h" +namespace Rodin::FormLanguage +{ + /** + * @ingroup TraitsSpecializations + */ + template + struct Traits> + { + using RangeType = Range; + }; +} + namespace Rodin::Variational { /**