Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
SuperellipseGenerator.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, are permitted provided that the following conditions are met:
14 * Redistributions of source code must retain the above copyright notice,
15 this list of conditions and the following disclaimer.
16 * Redistributions in binary form must reproduce the above copyright notice,
17 this list of conditions and the following disclaimer in the documentation
18 and/or other materials provided with the distribution.
19 * Neither the name of the University of Oxford nor the names of its
20 contributors may be used to endorse or promote products derived from this
21 software without specific prior written permission.
22
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36#include "SuperellipseGenerator.hpp"
37
39 double ellipseExponent,
40 double width,
41 double height,
42 double botLeftX,
43 double botLeftY)
44 : mTargetNodeSpacing(DOUBLE_UNSET),
45 mHeightOfTopSurface(0.3 * height)
46{
47 // Validate input
48 assert(numPoints > 1);
49 assert(ellipseExponent > 0.0);
50 assert(width > 0.0);
51 assert(height > 0.0);
52
53 // Set up member variables
54 mPoints.resize(numPoints);
55
56 /*
57 * Run a high-density pass around the parametric curve first
58 */
59 unsigned dense_pts = 50000;
60
61 // Vectors to store information from the loop
62 std::vector<double> cumulative_arc_length(dense_pts);
63 std::vector<c_vector<double, 2> > dense_locations(dense_pts);
64
65 // Helper variables for the loop
66 double dense_spacing = (2.0 * M_PI) / (double)dense_pts;
67 double cumulative_dist = 0.0;
68 double temp_sin, temp_cos;
69
70 // Fill in first location by hand
71 cumulative_arc_length[0] = 0.0;
72 dense_locations[0][0] = width * 0.5;
73 dense_locations[0][1] = 0.0;
74
75 // Fill in all other locations
76 for (unsigned point = 1; point < dense_pts; point++)
77 {
78 // Update cumulative distance
79 cumulative_dist += dense_spacing;
80
81 // Get the current sin and cos values
82 temp_cos = cos(cumulative_dist);
83 temp_sin = sin(cumulative_dist);
84
85 // x and y are powers of cos and sin, with the sign of sin and cos
86 dense_locations[point][0] = copysign(pow(fabs(temp_cos), ellipseExponent), temp_cos) * width * 0.5;
87 dense_locations[point][1] = copysign(pow(fabs(temp_sin), ellipseExponent), temp_sin) * height * 0.5;
88
89 // Check that the new point created lies within [-width/2, width/2] x [-height/2, height/2]
90 assert(fabs(dense_locations[point][0]) < width * 0.5 + 1e-10);
91 assert(fabs(dense_locations[point][1]) < height * 0.5 + 1e-10);
92
93 // Fill in arc length property
94 cumulative_arc_length[point] = cumulative_arc_length[point - 1] + norm_2(dense_locations[point] - dense_locations[point - 1]);
95 }
96
97 double total_arc_length = cumulative_arc_length[dense_pts - 1] + norm_2(dense_locations[dense_pts - 1] - dense_locations[0]);
98
99 // Since our perimeter is periodic, we get back to the start
100 cumulative_arc_length.push_back(total_arc_length);
101 dense_locations.push_back(dense_locations[0]);
102
103 /*
104 * Decide on the best place to put the actual boundary points
105 */
106
107 // Helper variables for loop
108 unsigned dense_it = 0;
109
110 mTargetNodeSpacing = total_arc_length / double(numPoints);
111
112 double target_arclength = 0.0;
113 double interpolant;
114
115 // Fill in first point by hand
116 mPoints[0] = dense_locations[0];
117
118 // Fill in all other locations
119 for (unsigned point = 1; point < numPoints; point++)
120 {
121 target_arclength = double(point) * mTargetNodeSpacing;
122
123 while (cumulative_arc_length[dense_it] < target_arclength && dense_it < dense_pts)
124 {
125 dense_it ++;
126 }
127
128 // Check we're one position either side of the target arclength
129 assert(cumulative_arc_length[dense_it - 1] < target_arclength);
130 assert(cumulative_arc_length[dense_it] >= target_arclength);
131
132 // Find the proportion between the two known locations we are
133 interpolant = (target_arclength - cumulative_arc_length[dense_it - 1]) / (cumulative_arc_length[dense_it] - cumulative_arc_length[dense_it - 1]);
134
135 // Calculate the location of the new point
136 mPoints[point] = (1.0 - interpolant) * dense_locations[dense_it - 1] + interpolant * dense_locations[dense_it];
137 }
138
139 /*
140 * If the exponent is at least 1, we return the default value which is the y-coord with largest value.
141 *
142 * If the exponent is less than 1, the top surface is defined as the height at which maximum curvature is attained.
143 * We loop through the points looking at the change in x and change in y from the previous point. The first
144 * difference in which the absolute change in x is greater than the absolute change in y will be the point with
145 * maximal curvature, for a superellipse.
146 */
147 if (ellipseExponent < 1.0)
148 {
149 double delta_x = fabs(mPoints[0][0] - mPoints[1][0]);
150 double delta_y = fabs(mPoints[0][1] - mPoints[1][1]);
151
152 // The following should always hold, for an exponent less than 1.0
153 if (delta_x > delta_y)
154 {
156 }
157
158 for (unsigned i = 1; i < numPoints; i++)
159 {
160 delta_x = fabs(mPoints[i][0] - mPoints[i - 1][0]);
161 delta_y = fabs(mPoints[i][1] - mPoints[i - 1][1]);
162
163 if (delta_x > delta_y)
164 {
165 mHeightOfTopSurface = 0.5 * (mPoints[i][1] + mPoints[i - 1][1]);
166
167 break;
168 }
169
170 // We should meet this condition before being half way through the list
171 if (i == unsigned(numPoints / 2.0))
172 {
174 }
175 }
176 }
177
178 /*
179 * Rescale all points to match parameters
180 */
181
182 // Move perimeter from [-width/2, width/2] x [-height/2, height/2] to correct location
183 c_vector<double, 2> offset;
184 offset[0] = width * 0.5 + botLeftX;
185 offset[1] = height * 0.5 + botLeftY;
186
187 // Reposition the height of the top surface
188 mHeightOfTopSurface += offset[1];
189
190 for (unsigned point = 0; point < numPoints; point++)
191 {
192 // Reposition all points
193 mPoints[point] += offset;
194
195 // Check we're in the right place
196 assert((mPoints[point][0] > botLeftX - 1e-10) && (mPoints[point][0] < botLeftX + width + 1e-10));
197 assert((mPoints[point][1] > botLeftY - 1e-10) && (mPoints[point][1] < botLeftY + height + 1e-10));
198 }
199}
200
205
210
215
216const std::vector<c_vector<double, 2> > SuperellipseGenerator::GetPointsAsVectors() const
217{
218 return mPoints;
219}
220
221const std::vector<ChastePoint<2> > SuperellipseGenerator::GetPointsAsChastePoints() const
222{
223 std::vector<ChastePoint<2> > chaste_points;
224
225 for (unsigned point = 0; point < mPoints.size(); point++)
226 {
227 chaste_points.push_back(ChastePoint<2>(mPoints[point]));
228 }
229
230 return chaste_points;
231}
const double DOUBLE_UNSET
Definition Exception.hpp:57
#define NEVER_REACHED
SuperellipseGenerator(unsigned numPoints=100, double ellipseExponent=1.0, double width=1.0, double height=1.0, double botLeftX=0.0, double botLeftY=0.0)
std::vector< c_vector< double, 2 > > mPoints
const std::vector< ChastePoint< 2 > > GetPointsAsChastePoints() const
const std::vector< c_vector< double, 2 > > GetPointsAsVectors() const