MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
material_metrics.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details
11
12#include "material_metrics.hpp"
13
14namespace mfem
15{
16
18{
19 std::vector<real_t> dist_vector;
20 dist_vector.resize(particle_positions_.size());
21 // 1. Compute the distance to each particle.
22 for (size_t i = 0; i < particle_positions_.size(); i++)
23 {
24 Vector y(3);
25 particle_orientations_[i].Mult(x, y);
26 dist_vector[i] = particle_positions_[i].DistanceTo(y);
27 }
28 // 2. Choose smallest number in the vector dist_vector.
29 real_t min_dist = *std::min_element(dist_vector.begin(), dist_vector.end());
30 return min_dist;
31}
32
33void ParticleTopology::Initialize(std::vector<real_t> &random_positions,
34 std::vector<real_t> &random_rotations)
35{
36 // 1. Initialize the particle positions.
37 particle_positions_.resize(number_of_particles_);
38 particle_orientations_.resize(number_of_particles_);
39 for (size_t i = 0; i < number_of_particles_; i++)
40 {
41 // 2.1 Read the positions.
42 size_t idx_pos = i * 3;
43 Vector particle_position({random_positions[idx_pos],
44 random_positions[idx_pos + 1],
45 random_positions[idx_pos + 2]});
46
47 // 2.2 Read the random rotations.
48 size_t idx_rot = i * 9;
49 DenseMatrix R(3, 3);
50 R(0, 0) = random_rotations[idx_rot + 0];
51 R(0, 1) = random_rotations[idx_rot + 1];
52 R(0, 2) = random_rotations[idx_rot + 2];
53 R(1, 0) = random_rotations[idx_rot + 3];
54 R(1, 1) = random_rotations[idx_rot + 4];
55 R(1, 2) = random_rotations[idx_rot + 5];
56 R(2, 0) = random_rotations[idx_rot + 6];
57 R(2, 1) = random_rotations[idx_rot + 7];
58 R(2, 2) = random_rotations[idx_rot + 8];
59
60 // 2.3 Fill the orientation vector.
61 DenseMatrix res(3, 3);
62 MultADBt(R, particle_shape_, R, res);
63 particle_orientations_[i] = res;
64
65 // 2.4 Scale position for distance metric
66 Vector scaled_position(3);
67 res.Mult(particle_position, scaled_position);
68 particle_positions_[i] = scaled_position;
69 }
70}
71
73{
74 // Implements formula used in [2, Example 5].
75 const real_t a = start_.DistanceTo(x);
76 const real_t b = end_.DistanceTo(x);
77 const real_t c = start_.DistanceTo(end_);
78 const real_t s1 = (pow(a, 2) + pow(b, 2)) / 2;
79 const real_t s2 = pow(c, 2) / 4;
80 const real_t s3 = pow((pow(a, 2) - pow(b, 2)) / (2 * c), 2);
81 return sqrt(abs(s1 - s2 - s3));
82}
83
85{
86 // Define the point in the vector which differentiates between the periodic
87 // points and the inner points.
88 constexpr size_t periodic_edges = 6;
89
90 // 1. Fill a vector with x and it's ghost points mimicking the periodicity
91 // of the topology.
92 std::vector<Vector> periodic_points;
93 CreatePeriodicPoints(x, periodic_points);
94 std::vector<real_t> dist_vector;
95
96 // 2. Compute the distance to each periodic points to the outer edges.
97 for (const auto &point : periodic_points)
98 {
99 for (size_t i = 0; i < periodic_edges; i++)
100 {
101 dist_vector.push_back(edges_[i].GetDistanceTo(point));
102 }
103 }
104
105 // 3. Add distance between x and the remaining inner edges
106 for (size_t i = periodic_edges; i < edges_.size(); i++)
107 {
108 dist_vector.push_back(edges_[i].GetDistanceTo(x));
109 }
110
111 // 3. Choose the smallest number in the vector dist_vector.
112 real_t min_dist = *std::min_element(dist_vector.begin(), dist_vector.end());
113 return min_dist;
114}
115
116void OctetTrussTopology::Initialize()
117{
118 // 1. Create the points defining the topology (begin and end points of the
119 // edges). Outer structure
120 real_t p1_data[3] = {0, 0, 0};
121 real_t p2_data[3] = {0, 1, 1};
122 real_t p3_data[3] = {1, 0, 1};
123 real_t p4_data[3] = {1, 1, 0};
124
125 Vector p1(p1_data, 3);
126 Vector p2(p2_data, 3);
127 Vector p3(p3_data, 3);
128 Vector p4(p4_data, 3);
129
130 points_.push_back(p1);
131 points_.push_back(p2);
132 points_.push_back(p3);
133 points_.push_back(p4);
134
135 // 2. Create the inner structure
136 real_t p5_data[3] = {0, 0.5, 0.5}; // left
137 real_t p6_data[3] = {1, 0.5, 0.5}; // right
138 real_t p7_data[3] = {0.5, 0, 0.5}; // bottom
139 real_t p8_data[3] = {0.5, 1, 0.5}; // top
140 real_t p9_data[3] = {0.5, 0.5, 0}; // front
141 real_t p10_data[3] = {0.5, 0.5, 1}; // back
142
143 Vector p5(p5_data, 3);
144 Vector p6(p6_data, 3);
145 Vector p7(p7_data, 3);
146 Vector p8(p8_data, 3);
147 Vector p9(p9_data, 3);
148 Vector p10(p10_data, 3);
149
150 points_.push_back(p5);
151 points_.push_back(p6);
152 points_.push_back(p7);
153 points_.push_back(p8);
154 points_.push_back(p9);
155 points_.push_back(p10);
156
157 // 3. Create the outer edges.
158 for (size_t i = 0; i < 4; i++)
159 {
160 for (size_t j = i + 1; j < 4; j++)
161 {
162 Edge edge(points_[i], points_[j]);
163 edges_.push_back(edge);
164 }
165 }
166
167 // 4. Create the inner edges from p5 and p6 to p7, p8, p9, p10; plus the
168 // latter four connected in a cycle.
169 Edge edge5(points_[4], points_[6]);
170 Edge edge6(points_[4], points_[7]);
171 Edge edge7(points_[4], points_[8]);
172 Edge edge8(points_[4], points_[9]);
173 Edge edge9(points_[5], points_[6]);
174 Edge edge10(points_[5], points_[7]);
175 Edge edge11(points_[5], points_[8]);
176 Edge edge12(points_[5], points_[9]);
177 Edge edge13(points_[6], points_[8]);
178 Edge edge14(points_[6], points_[9]);
179 Edge edge15(points_[7], points_[8]);
180 Edge edge16(points_[7], points_[9]);
181
182 // Push into edges vector
183 edges_.push_back(edge5);
184 edges_.push_back(edge6);
185 edges_.push_back(edge7);
186 edges_.push_back(edge8);
187 edges_.push_back(edge9);
188 edges_.push_back(edge10);
189 edges_.push_back(edge11);
190 edges_.push_back(edge12);
191 edges_.push_back(edge13);
192 edges_.push_back(edge14);
193 edges_.push_back(edge15);
194 edges_.push_back(edge16);
195}
196
197void OctetTrussTopology::CreatePeriodicPoints(
198 const Vector &x, std::vector<Vector> &periodic_points) const
199{
200 Vector xx(x);
201 // Compute the displaced ghost points. Computation assumes domain [0,1]^3.
202 real_t d_x[3] = {1, 0, 0};
203 real_t d_y[3] = {0, 1, 0};
204 real_t d_z[3] = {0, 0, 1};
205
206 Vector dispcement_x(d_x, 3);
207 Vector dispcement_y(d_y, 3);
208 Vector dispcement_z(d_z, 3);
209
210 Vector x_shifted_x_pos = x;
211 x_shifted_x_pos += dispcement_x;
212 Vector x_shifted_x_neg = x;
213 x_shifted_x_neg -= dispcement_x;
214 Vector x_shifted_y_pos = x;
215 x_shifted_y_pos += dispcement_y;
216 Vector x_shifted_y_neg = x;
217 x_shifted_y_neg -= dispcement_y;
218 Vector x_shifted_z_pos = x;
219 x_shifted_z_pos += dispcement_z;
220 Vector x_shifted_z_neg = x;
221 x_shifted_z_neg -= dispcement_z;
222 // Fill the vector with all relevant points
223 periodic_points.push_back(xx);
224 periodic_points.push_back(x_shifted_x_pos);
225 periodic_points.push_back(x_shifted_x_neg);
226 periodic_points.push_back(x_shifted_y_pos);
227 periodic_points.push_back(x_shifted_y_neg);
228 periodic_points.push_back(x_shifted_z_pos);
229 periodic_points.push_back(x_shifted_z_neg);
230}
231
232} // namespace mfem
real_t GetDistanceTo(const Vector &x) const
Compute the distance between a point and the edge.
real_t ComputeMetric(const Vector &x) final
Compute the metric rho describing the material topology.
real_t ComputeMetric(const Vector &x) final
Vector data type.
Definition vector.hpp:80
real_t DistanceTo(const real_t *p) const
Compute the Euclidean distance to another vector.
Definition vector.hpp:690
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
void MultADBt(const DenseMatrix &A, const Vector &D, const DenseMatrix &B, DenseMatrix &ADBt)
ADBt = A D B^t, where D is diagonal.
float real_t
Definition config.hpp:43