[cig-commits] [commit] master: astyle cleanup (a11da16)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Oct 8 17:05:08 PDT 2014
Repository : https://github.com/geodynamics/vc
On branch : master
Link : https://github.com/geodynamics/vc/compare/23464fca3efa2b6ad7ee0ce8f60c225b18b49741...e4325192ad1118379f46ba66899cb98143d09e04
>---------------------------------------------------------------
commit a11da169e1b1b45b62a1d09c6831d342069c01ff
Author: Eric Heien <emheien at ucdavis.edu>
Date: Wed Sep 17 22:27:48 2014 -0700
astyle cleanup
>---------------------------------------------------------------
a11da169e1b1b45b62a1d09c6831d342069c01ff
quakelib/src/QuakeLibIO.cpp | 203 +++++++++++++++++++----------------
src/core/VCEvent.h | 0
src/core/VCSimDataEvents.h | 0
src/io/EventOutput.cpp | 1 +
src/mesher.cpp | 0
src/simulation/BASSAftershocks.cpp | 2 +-
src/simulation/BASSAftershocks.h | 0
src/simulation/UpdateBlockStress.cpp | 0
8 files changed, 112 insertions(+), 94 deletions(-)
diff --git a/quakelib/src/QuakeLibIO.cpp b/quakelib/src/QuakeLibIO.cpp
index 09c3cec..6641216 100644
--- a/quakelib/src/QuakeLibIO.cpp
+++ b/quakelib/src/QuakeLibIO.cpp
@@ -334,102 +334,112 @@ void quakelib::FaultTracePoint::write_ascii(std::ostream &out_stream) const {
representing how long along the trace a point is. The spline is a linear interpolation between successive points.
*/
class TraceSpline {
-private:
- //! The points comprising the spline
- std::vector<quakelib::Vec<3> > _pts;
- //! Total length of the distances between successive points
- double _spline_len;
- //! Individual lengths between successive points
- //! _point_dists[i] is the distance between _pts[i] and _pts[i+1]
- std::vector<double> _point_dists;
-
-public:
- TraceSpline(void) : _spline_len(0) {};
-
- // Add another point to the spline
- void add_point(const quakelib::Vec<3> &new_pt) {
- double add_dist = (_pts.size() > 0 ? _pts.back().dist(new_pt) : 0);
- _spline_len += add_dist;
- if (_pts.size() > 0) _point_dists.push_back(add_dist);
- _pts.push_back(new_pt);
- }
-
- // Return the element index and inner t corresponding to parameter t
- void get_element(const double t, unsigned int &index, double &inner_t) {
- double spline_dist = t * _spline_len;
-
- // Ensure t is in [0,1]
- assertThrow(t >= 0 && t <= 1, std::out_of_range("TraceSpline::interpolate"));
-
- // Go through each point
- for (unsigned int i=0;i<_pts.size()-1;++i) {
- // If we're between the points for this t, interpolate and return
- if (spline_dist <= _point_dists[i]) {
- index = i;
- if (_point_dists[i] != 0) inner_t = spline_dist / _point_dists[i];
- else inner_t = 0;
- return;
- }
- spline_dist -= _point_dists[i];
+ private:
+ //! The points comprising the spline
+ std::vector<quakelib::Vec<3> > _pts;
+ //! Total length of the distances between successive points
+ double _spline_len;
+ //! Individual lengths between successive points
+ //! _point_dists[i] is the distance between _pts[i] and _pts[i+1]
+ std::vector<double> _point_dists;
+
+ public:
+ TraceSpline(void) : _spline_len(0) {};
+
+ // Add another point to the spline
+ void add_point(const quakelib::Vec<3> &new_pt) {
+ double add_dist = (_pts.size() > 0 ? _pts.back().dist(new_pt) : 0);
+ _spline_len += add_dist;
+
+ if (_pts.size() > 0) _point_dists.push_back(add_dist);
+
+ _pts.push_back(new_pt);
}
- // If we reach the end, we return the final point and inner_t = 0
- index = _pts.size()-1;
- inner_t = 0;
- }
- // Return the point on this spline at t
- quakelib::Vec<3> interpolate(const double t) {
- unsigned int ind;
- double inner_t;
+ // Return the element index and inner t corresponding to parameter t
+ void get_element(const double t, unsigned int &index, double &inner_t) {
+ double spline_dist = t * _spline_len;
- assertThrow(t >= 0 && t <= 1, std::out_of_range("TraceSpline::interpolate"));
- get_element(t, ind, inner_t);
- if (ind == _pts.size() - 1) return _pts.back();
- else return _pts.at(ind) * (1-inner_t) + _pts.at(ind+1) * inner_t;
- }
+ // Ensure t is in [0,1]
+ assertThrow(t >= 0 && t <= 1, std::out_of_range("TraceSpline::interpolate"));
- // Given a starting point t (in [0,1]) and an element size,
- // returns the next t which is linear elem_size away (not distance along the spline)
- double advance_element(const double t, const double elem_size) {
- unsigned int ind;
- double inner_t;
+ // Go through each point
+ for (unsigned int i=0; i<_pts.size()-1; ++i) {
+ // If we're between the points for this t, interpolate and return
+ if (spline_dist <= _point_dists[i]) {
+ index = i;
- // Find the starting point on the spline at t
- quakelib::Vec<3> start_pt = interpolate(t);
+ if (_point_dists[i] != 0) inner_t = spline_dist / _point_dists[i];
+ else inner_t = 0;
- // Get the element that the starting point is associated iwth
- get_element(t, ind, inner_t);
+ return;
+ }
- // Keep going until we find a trace point which would create
- // an element greater than our target size
- double cur_dist = t * _spline_len;
- while (ind+1 < _pts.size() && start_pt.dist(_pts.at(ind+1)) < elem_size) {
- ind++;
- if (ind < _point_dists.size()) cur_dist += (1-inner_t) * _point_dists.at(ind);
- else cur_dist += elem_size;
+ spline_dist -= _point_dists[i];
+ }
+
+ // If we reach the end, we return the final point and inner_t = 0
+ index = _pts.size()-1;
inner_t = 0;
}
- // If we're past the end of the trace, return our best guess
- // for t based on the size of the last segment
- // This is needed to adjust the element size during meshing
- if (ind+1 == _pts.size()) return cur_dist/_spline_len;
-
- // Now we know the points between which our element must exist (ind and ind+1)
- double x_0 = _pts.at(ind)[0], x_1 = _pts.at(ind+1)[0], x_s = start_pt[0];
- double y_0 = _pts.at(ind)[1], y_1 = _pts.at(ind+1)[1], y_s = start_pt[1], l = elem_size;
-
- // Calculate the inner t between these two points
- double next_t = fabs((sqrt(pow(-2*x_0*x_0+2*x_0*x_1+2*x_0*x_s-2*x_1*x_s-2*y_0*y_0+2*y_0*y_1+2*y_0*y_s-2*y_1*y_s, 2)-
- 4*(x_0*x_0-2*x_0*x_1+x_1*x_1+y_0*y_0-2*y_0*y_1+y_1*y_1)*
- (x_0*x_0-2*x_0*x_s+x_s*x_s+y_0*y_0-2*y_0*y_s+y_s*y_s-l*l))+
- 2*x_0*x_0-2*x_0*x_1-2*x_0*x_s+2*x_1*x_s+2*y_0*y_0-2*y_0*y_1-2*y_0*y_s+2*y_1*y_s)/
- (2*(x_0*x_0-2*x_0*x_1+x_1*x_1+y_0*y_0-2*y_0*y_1+y_1*y_1)));
-
- // Given this point, recalculate t and return
- cur_dist += (next_t-inner_t) * _point_dists[ind];
- return cur_dist/_spline_len;
- };
+ // Return the point on this spline at t
+ quakelib::Vec<3> interpolate(const double t) {
+ unsigned int ind;
+ double inner_t;
+
+ assertThrow(t >= 0 && t <= 1, std::out_of_range("TraceSpline::interpolate"));
+ get_element(t, ind, inner_t);
+
+ if (ind == _pts.size() - 1) return _pts.back();
+ else return _pts.at(ind) * (1-inner_t) + _pts.at(ind+1) * inner_t;
+ }
+
+ // Given a starting point t (in [0,1]) and an element size,
+ // returns the next t which is linear elem_size away (not distance along the spline)
+ double advance_element(const double t, const double elem_size) {
+ unsigned int ind;
+ double inner_t;
+
+ // Find the starting point on the spline at t
+ quakelib::Vec<3> start_pt = interpolate(t);
+
+ // Get the element that the starting point is associated iwth
+ get_element(t, ind, inner_t);
+
+ // Keep going until we find a trace point which would create
+ // an element greater than our target size
+ double cur_dist = t * _spline_len;
+
+ while (ind+1 < _pts.size() && start_pt.dist(_pts.at(ind+1)) < elem_size) {
+ ind++;
+
+ if (ind < _point_dists.size()) cur_dist += (1-inner_t) * _point_dists.at(ind);
+ else cur_dist += elem_size;
+
+ inner_t = 0;
+ }
+
+ // If we're past the end of the trace, return our best guess
+ // for t based on the size of the last segment
+ // This is needed to adjust the element size during meshing
+ if (ind+1 == _pts.size()) return cur_dist/_spline_len;
+
+ // Now we know the points between which our element must exist (ind and ind+1)
+ double x_0 = _pts.at(ind)[0], x_1 = _pts.at(ind+1)[0], x_s = start_pt[0];
+ double y_0 = _pts.at(ind)[1], y_1 = _pts.at(ind+1)[1], y_s = start_pt[1], l = elem_size;
+
+ // Calculate the inner t between these two points
+ double next_t = fabs((sqrt(pow(-2*x_0*x_0+2*x_0*x_1+2*x_0*x_s-2*x_1*x_s-2*y_0*y_0+2*y_0*y_1+2*y_0*y_s-2*y_1*y_s, 2)-
+ 4*(x_0*x_0-2*x_0*x_1+x_1*x_1+y_0*y_0-2*y_0*y_1+y_1*y_1)*
+ (x_0*x_0-2*x_0*x_s+x_s*x_s+y_0*y_0-2*y_0*y_s+y_s*y_s-l*l))+
+ 2*x_0*x_0-2*x_0*x_1-2*x_0*x_s+2*x_1*x_s+2*y_0*y_0-2*y_0*y_1-2*y_0*y_s+2*y_1*y_s)/
+ (2*(x_0*x_0-2*x_0*x_1+x_1*x_1+y_0*y_0-2*y_0*y_1+y_1*y_1)));
+
+ // Given this point, recalculate t and return
+ cur_dist += (next_t-inner_t) * _point_dists[ind];
+ return cur_dist/_spline_len;
+ };
};
void quakelib::ModelWorld::create_section(std::vector<unsigned int> &unused_trace_segments, const std::vector<FaultTracePoint> &trace, const LatLonDepth &base_coord, const UIndex &fault_id, const float &element_size, const std::string §ion_name, const std::string &taper_method) {
@@ -448,6 +458,7 @@ void quakelib::ModelWorld::create_section(std::vector<unsigned int> &unused_trac
if (element_size <= 0) return;
num_trace_pts = trace.size();
+
if (num_trace_pts == 0 || num_trace_pts == 1) return;
ModelSection §ion = new_section();
@@ -455,7 +466,7 @@ void quakelib::ModelWorld::create_section(std::vector<unsigned int> &unused_trac
section.set_fault_id(fault_id);
// Create a spline with the trace points
- for (i=0;i<num_trace_pts;++i) {
+ for (i=0; i<num_trace_pts; ++i) {
Vec<3> pt = conv.convert2xyz(trace.at(i).pos());
spline.add_point(pt);
unused_trace_pts.insert(i);
@@ -468,12 +479,15 @@ void quakelib::ModelWorld::create_section(std::vector<unsigned int> &unused_trac
double cur_elem_size_guess = element_size;
double overshoot = DBL_MAX, undershoot = DBL_MAX, t_per_meter=1;
unsigned int num_iters = 0;
+
while ((fmin(overshoot, undershoot)/t_per_meter)/element_size > 1e-3 && num_iters < 10) {
double cur_t, next_t, sum_t=0;
elem_count = 0;
cur_t = 0;
+
while (cur_t < 1) {
next_t = spline.advance_element(cur_t, cur_elem_size_guess);
+
if (next_t < 1) {
sum_t += next_t-cur_t;
elem_count++;
@@ -482,13 +496,16 @@ void quakelib::ModelWorld::create_section(std::vector<unsigned int> &unused_trac
break;
}
}
+
t_per_meter = (sum_t/elem_count)/cur_elem_size_guess;
overshoot = fabs(1-next_t)/elem_count;
undershoot = fabs(1-cur_t)/elem_count;
+
// If we're undershooting by less than we're overshooting, the best solution is to increase block size
if (undershoot < overshoot) cur_elem_size_guess += undershoot/t_per_meter;
// Otherwise we decrease block size
else cur_elem_size_guess -= overshoot/t_per_meter;
+
num_iters++;
}
@@ -505,7 +522,7 @@ void quakelib::ModelWorld::create_section(std::vector<unsigned int> &unused_trac
cur_elem_ind = 0;
unused_trace_pts.erase(0);
- for (i=0;i<elem_count;++i) {
+ for (i=0; i<elem_count; ++i) {
// Get the next t value along the trace by advancing the element size
next_t = spline.advance_element(cur_t, horiz_elem_size);
@@ -534,12 +551,12 @@ void quakelib::ModelWorld::create_section(std::vector<unsigned int> &unused_trac
if (num_vert_elems == 0) std::cerr << "WARNING: Depth is smaller than element size in trace segment "
<< next_elem_ind << ". Element size may be too big." << std::endl;
- elem_slip_rate = conv.cm_per_yr2m_per_sec(inner_t*trace.at(next_elem_ind).slip_rate()+(1.0-inner_t)*trace.at(cur_elem_ind).slip_rate());
- elem_aseismic = inner_t*trace.at(next_elem_ind).aseismic()+(1.0-inner_t)*trace.at(cur_elem_ind).aseismic();
- elem_dip = conv.deg2rad(inner_t*trace.at(next_elem_ind).dip()+(1.0-inner_t)*trace.at(cur_elem_ind).dip());
- elem_rake = conv.deg2rad(inner_t*trace.at(next_elem_ind).rake()+(1.0-inner_t)*trace.at(cur_elem_ind).rake());
+ elem_slip_rate = conv.cm_per_yr2m_per_sec(inner_t *trace.at(next_elem_ind).slip_rate()+(1.0-inner_t)*trace.at(cur_elem_ind).slip_rate());
+ elem_aseismic = inner_t *trace.at(next_elem_ind).aseismic()+(1.0-inner_t)*trace.at(cur_elem_ind).aseismic();
+ elem_dip = conv.deg2rad(inner_t *trace.at(next_elem_ind).dip()+(1.0-inner_t)*trace.at(cur_elem_ind).dip());
+ elem_rake = conv.deg2rad(inner_t *trace.at(next_elem_ind).rake()+(1.0-inner_t)*trace.at(cur_elem_ind).rake());
elem_lame_mu = inner_t *trace.at(next_elem_ind).lame_mu()+(1.0-inner_t)*trace.at(cur_elem_ind).lame_mu();
- elem_lame_lambda = inner_t*trace.at(next_elem_ind).lame_lambda()+(1.0-inner_t)*trace.at(cur_elem_ind).lame_lambda();
+ elem_lame_lambda = inner_t *trace.at(next_elem_ind).lame_lambda()+(1.0-inner_t)*trace.at(cur_elem_ind).lame_lambda();
// Set up the vertical step to go along the dip
vert_step = element_step_vec.rotate_around_axis(Vec<3>(0,0,-1), M_PI/2);
diff --git a/src/io/EventOutput.cpp b/src/io/EventOutput.cpp
index 1507e9f..16b3285 100644
--- a/src/io/EventOutput.cpp
+++ b/src/io/EventOutput.cpp
@@ -91,6 +91,7 @@ void EventOutput::init(SimFramework *_sim) {
std::cerr << "ERROR: Could not open output file " << sim->getEventOutfile() << std::endl;
exit(-1);
}
+
writeEventFileHeader();
writeSweepFileHeader();
} else {
diff --git a/src/simulation/BASSAftershocks.cpp b/src/simulation/BASSAftershocks.cpp
index a5db798..b36173d 100644
--- a/src/simulation/BASSAftershocks.cpp
+++ b/src/simulation/BASSAftershocks.cpp
@@ -82,7 +82,7 @@ SimRequest BASSAftershocks::run(SimFramework *_sim) {
// Finally we sort the events into time order since the aftershocks in the
// BASS model are usually not generated in ordered time sequence.
- for (it=events_to_process.begin();it!=events_to_process.end();++it) {
+ for (it=events_to_process.begin(); it!=events_to_process.end(); ++it) {
sim->addAftershock(*it);
}
More information about the CIG-COMMITS
mailing list