|
| 1 | +// SPDX-FileCopyrightText: 2025 VTT Technical Research Centre of Finland Ltd |
| 2 | +// SPDX-License-Identifier: AGPL-3.0-or-later |
| 3 | + |
| 4 | +/** |
| 5 | + * @file 14_custom_field_initializer.cpp |
| 6 | + * @brief Example: Custom field initialization patterns |
| 7 | + * |
| 8 | + * @details |
| 9 | + * This example demonstrates how to create custom field initialization patterns |
| 10 | + * using simple structs and evaluation functions. We show three physical patterns: |
| 11 | + * |
| 12 | + * 1. **Lamb-Oseen Vortex** - Rotating fluid vortex with viscous core |
| 13 | + * 2. **Gaussian Bump** - Localized concentration or temperature field |
| 14 | + * 3. **Checkerboard** - Periodic alternating pattern |
| 15 | + * |
| 16 | + * ## Key Concept: Struct-Based Patterns |
| 17 | + * |
| 18 | + * Instead of inheritance hierarchies, we use simple structs with public data. |
| 19 | + * This follows OpenPFC's "structs + free functions" philosophy. |
| 20 | + * |
| 21 | + * ## Integration with OpenPFC |
| 22 | + * |
| 23 | + * For full integration with DiscreteField/World APIs using ADL, see: |
| 24 | + * - docs/extending_openpfc/adl_extension_patterns.md |
| 25 | + * - examples/17_custom_coordinate_system.cpp |
| 26 | + * |
| 27 | + * This simplified example focuses on the pattern concept itself. |
| 28 | + */ |
| 29 | + |
| 30 | +#include <openpfc/openpfc.hpp> |
| 31 | +#include <cmath> |
| 32 | +#include <iostream> |
| 33 | + |
| 34 | +using namespace pfc; |
| 35 | + |
| 36 | +// Use M_PI constant |
| 37 | +#ifndef M_PI |
| 38 | +constexpr double M_PI = 3.14159265358979323846; |
| 39 | +#endif |
| 40 | + |
| 41 | +// ============================================================================ |
| 42 | +// Part 1: Define Custom Pattern Types |
| 43 | +// ============================================================================ |
| 44 | + |
| 45 | +/** |
| 46 | + * Custom namespace for user-defined patterns. |
| 47 | + */ |
| 48 | +namespace my_project { |
| 49 | + |
| 50 | +/** |
| 51 | + * @brief Lamb-Oseen vortex pattern |
| 52 | + * |
| 53 | + * Models a rotating vortex with viscous core, common in fluid dynamics. |
| 54 | + * The tangential velocity follows: v_θ(r) = (Γ/2πr)[1 - exp(-r²/r_c²)] |
| 55 | + * |
| 56 | + * @see https://en.wikipedia.org/wiki/Lamb%E2%80%93Oseen_vortex |
| 57 | + */ |
| 58 | +struct VortexPattern { |
| 59 | + Real3 m_center; ///< Vortex center (x, y, z) |
| 60 | + double m_strength; ///< Circulation Γ |
| 61 | + double m_core_radius; ///< Core radius r_c |
| 62 | + |
| 63 | + VortexPattern(Real3 center, double strength, double core_radius) |
| 64 | + : m_center(center), m_strength(strength), m_core_radius(core_radius) {} |
| 65 | +}; |
| 66 | + |
| 67 | +/** |
| 68 | + * @brief 3D Gaussian bump |
| 69 | + * |
| 70 | + * Models localized concentration, temperature, or density field. |
| 71 | + * φ(r) = A * exp(-r²/(2σ²)) |
| 72 | + */ |
| 73 | +struct GaussianBump { |
| 74 | + Real3 m_center; ///< Peak location |
| 75 | + double m_amplitude; ///< Peak height A |
| 76 | + double m_width; ///< Standard deviation σ |
| 77 | + |
| 78 | + GaussianBump(Real3 center, double amplitude, double width) |
| 79 | + : m_center(center), m_amplitude(amplitude), m_width(width) {} |
| 80 | +}; |
| 81 | + |
| 82 | +/** |
| 83 | + * @brief 3D checkerboard pattern |
| 84 | + * |
| 85 | + * Periodic alternating values, useful for testing and validation. |
| 86 | + */ |
| 87 | +struct CheckerboardPattern { |
| 88 | + double m_value_high; ///< Value in "white" cells |
| 89 | + double m_value_low; ///< Value in "black" cells |
| 90 | + Real3 m_period; ///< Period in each direction |
| 91 | + |
| 92 | + CheckerboardPattern(double high, double low, Real3 period) |
| 93 | + : m_value_high(high), m_value_low(low), m_period(period) {} |
| 94 | +}; |
| 95 | + |
| 96 | +} // namespace my_project |
| 97 | + |
| 98 | +// ============================================================================ |
| 99 | +// Part 2: Evaluation Functions |
| 100 | +// ============================================================================ |
| 101 | + |
| 102 | +/** |
| 103 | + * @brief Evaluate vortex pattern at given position |
| 104 | + * @param pattern The vortex configuration |
| 105 | + * @param pos Physical position to evaluate at |
| 106 | + * @return Tangential velocity at position |
| 107 | + */ |
| 108 | +double evaluate_vortex(const my_project::VortexPattern& pattern, const Real3& pos) { |
| 109 | + // Distance from vortex center in x-y plane |
| 110 | + double dx = pos[0] - pattern.m_center[0]; |
| 111 | + double dy = pos[1] - pattern.m_center[1]; |
| 112 | + double r = std::sqrt(dx * dx + dy * dy); |
| 113 | + |
| 114 | + // Lamb-Oseen vortex profile |
| 115 | + double r_c_sq = pattern.m_core_radius * pattern.m_core_radius; |
| 116 | + double value = 0.0; |
| 117 | + |
| 118 | + if (r > 1e-10) { // Avoid division by zero at center |
| 119 | + value = (pattern.m_strength / (2.0 * M_PI * r)) * |
| 120 | + (1.0 - std::exp(-r * r / r_c_sq)); |
| 121 | + } |
| 122 | + |
| 123 | + return value; |
| 124 | +} |
| 125 | + |
| 126 | +/** |
| 127 | + * @brief Evaluate Gaussian bump at given position |
| 128 | + * @param pattern The Gaussian configuration |
| 129 | + * @param pos Physical position to evaluate at |
| 130 | + * @return Field value at position |
| 131 | + */ |
| 132 | +double evaluate_gaussian(const my_project::GaussianBump& pattern, const Real3& pos) { |
| 133 | + // Distance from center |
| 134 | + double dx = pos[0] - pattern.m_center[0]; |
| 135 | + double dy = pos[1] - pattern.m_center[1]; |
| 136 | + double dz = pos[2] - pattern.m_center[2]; |
| 137 | + double dist_sq = dx*dx + dy*dy + dz*dz; |
| 138 | + |
| 139 | + // Gaussian: φ = A * exp(-dist² / (2σ²)) |
| 140 | + double sigma_sq = pattern.m_width * pattern.m_width; |
| 141 | + return pattern.m_amplitude * std::exp(-dist_sq / (2.0 * sigma_sq)); |
| 142 | +} |
| 143 | + |
| 144 | +/** |
| 145 | + * @brief Evaluate checkerboard at given position |
| 146 | + * @param pattern The checkerboard configuration |
| 147 | + * @param pos Physical position to evaluate at |
| 148 | + * @return High or low value depending on position |
| 149 | + */ |
| 150 | +double evaluate_checkerboard(const my_project::CheckerboardPattern& pattern, const Real3& pos) { |
| 151 | + // Determine which cell of the checkerboard |
| 152 | + int cell_i = static_cast<int>(std::floor(pos[0] / pattern.m_period[0])); |
| 153 | + int cell_j = static_cast<int>(std::floor(pos[1] / pattern.m_period[1])); |
| 154 | + int cell_k = static_cast<int>(std::floor(pos[2] / pattern.m_period[2])); |
| 155 | + |
| 156 | + // Checkerboard: alternate based on sum of cell indices |
| 157 | + int sum = cell_i + cell_j + cell_k; |
| 158 | + return (sum % 2 == 0) ? pattern.m_value_high : pattern.m_value_low; |
| 159 | +} |
| 160 | + |
| 161 | +// ============================================================================ |
| 162 | +// Part 3: Usage Examples |
| 163 | +// ============================================================================ |
| 164 | + |
| 165 | +void example_vortex_pattern() { |
| 166 | + std::cout << "=== Example 1: Vortex Pattern ===\n\n"; |
| 167 | + |
| 168 | + // Define vortex at origin with strength 5.0 and core radius 2.0 |
| 169 | + my_project::VortexPattern vortex({0.0, 0.0, 0.0}, 5.0, 2.0); |
| 170 | + |
| 171 | + // Evaluate at several radial distances |
| 172 | + std::cout << "Vortex tangential velocity profile:\n"; |
| 173 | + for (double r = 0.0; r <= 10.0; r += 2.0) { |
| 174 | + Real3 pos{r, 0.0, 0.0}; // Along x-axis |
| 175 | + double velocity = evaluate_vortex(vortex, pos); |
| 176 | + std::cout << " r = " << r << " : v_θ = " << velocity << "\n"; |
| 177 | + } |
| 178 | + std::cout << "\n"; |
| 179 | +} |
| 180 | + |
| 181 | +void example_gaussian_bump() { |
| 182 | + std::cout << "=== Example 2: Gaussian Bump ===\n\n"; |
| 183 | + |
| 184 | + // Create Gaussian peak at origin with amplitude 1.0, width 1.5 |
| 185 | + my_project::GaussianBump bump({0.0, 0.0, 0.0}, 1.0, 1.5); |
| 186 | + |
| 187 | + // Evaluate along a line |
| 188 | + std::cout << "Gaussian profile:\n"; |
| 189 | + for (double x = 0.0; x <= 5.0; x += 1.0) { |
| 190 | + Real3 pos{x, 0.0, 0.0}; |
| 191 | + double value = evaluate_gaussian(bump, pos); |
| 192 | + std::cout << " x = " << x << " : φ = " << value << "\n"; |
| 193 | + } |
| 194 | + std::cout << "\n"; |
| 195 | +} |
| 196 | + |
| 197 | +void example_checkerboard() { |
| 198 | + std::cout << "=== Example 3: Checkerboard Pattern ===\n\n"; |
| 199 | + |
| 200 | + // Create checkerboard with period 2.0 in each direction |
| 201 | + my_project::CheckerboardPattern checker(1.0, -1.0, {2.0, 2.0, 2.0}); |
| 202 | + |
| 203 | + // Sample a 2D slice (z=0 plane) |
| 204 | + std::cout << "Checkerboard pattern (z=0 plane):\n"; |
| 205 | + for (int j = 0; j < 4; ++j) { |
| 206 | + std::cout << " "; |
| 207 | + for (int i = 0; i < 4; ++i) { |
| 208 | + Real3 pos{i * 1.0, j * 1.0, 0.0}; |
| 209 | + double value = evaluate_checkerboard(checker, pos); |
| 210 | + std::cout << (value > 0 ? "+" : "-") << " "; |
| 211 | + } |
| 212 | + std::cout << "\n"; |
| 213 | + } |
| 214 | + std::cout << "\n"; |
| 215 | +} |
| 216 | + |
| 217 | +int main() { |
| 218 | + std::cout << "\n"; |
| 219 | + std::cout << "╔════════════════════════════════════════════════════════╗\n"; |
| 220 | + std::cout << "║ OpenPFC: Custom Field Initialization Patterns Example ║\n"; |
| 221 | + std::cout << "╚════════════════════════════════════════════════════════╝\n"; |
| 222 | + std::cout << "\n"; |
| 223 | + |
| 224 | + example_vortex_pattern(); |
| 225 | + example_gaussian_bump(); |
| 226 | + example_checkerboard(); |
| 227 | + |
| 228 | + std::cout << "✅ All examples completed successfully!\n"; |
| 229 | + std::cout << "\n"; |
| 230 | + std::cout << "📖 For full integration with DiscreteField/World using ADL, see:\n"; |
| 231 | + std::cout << " - docs/extending_openpfc/adl_extension_patterns.md\n"; |
| 232 | + std::cout << " - examples/17_custom_coordinate_system.cpp\n"; |
| 233 | + std::cout << "\n"; |
| 234 | + |
| 235 | + return 0; |
| 236 | +} |
0 commit comments