svMultiPhysics
Loading...
Searching...
No Matches
Tensor4.h
1/* Copyright (c) Stanford University, The Regents of the University of California, and others.
2 *
3 * All Rights Reserved.
4 *
5 * See Copyright-SimVascular.txt for additional details.
6 *
7 * Permission is hereby granted, free of charge, to any person obtaining
8 * a copy of this software and associated documentation files (the
9 * "Software"), to deal in the Software without restriction, including
10 * without limitation the rights to use, copy, modify, merge, publish,
11 * distribute, sublicense, and/or sell copies of the Software, and to
12 * permit persons to whom the Software is furnished to do so, subject
13 * to the following conditions:
14 *
15 * The above copyright notice and this permission notice shall be included
16 * in all copies or substantial portions of the Software.
17 *
18 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
19 * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
20 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
21 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
22 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
23 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
24 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
25 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
26 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
27 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
28 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 */
30
31#ifndef TENSOR4_H
32#define TENSOR4_H
33
34#include <cstring>
35#include <iostream>
36
37#ifdef ENABLE_ARRAY_INDEX_CHECKING
38#define Tensor4_check_enabled
39#endif
40
41/// @brief The Tensor4 template class implements a simple interface to 4th order tensors.
42//
43template<typename T>
44class Tensor4
45{
46 static bool show_index_check_message;
47
48 public:
49 std::string name_ = "";
50 int ni_ = 0;
51 int nj_ = 0;
52 int nk_ = 0;
53 int nl_ = 0;
54 int p1_ = 0;
55 int p2_ = 0;
56 int size_ = 0;
57 T *data_ = nullptr;
58
59 Tensor4()
60 {
61 std::string name_ = "";
62 int ni_ = 0;
63 int nj_ = 0;
64 int nk_ = 0;
65 int nl_ = 0;
66 int p1_ = 0;
67 int p2_ = 0;
68 int size_ = 0;
69 T *data_ = nullptr;
70 };
71
72 Tensor4(const int num_i, const int num_j, const int num_k, const int num_l)
73 {
74 allocate(num_i, num_j, num_k, num_l);
75 }
76
77 ~Tensor4()
78 {
79 //std::cout << "- - - - - Tensor4 dtor - - - - - " << std::endl;
80 if (data_ != nullptr) {
81 //std::cout << "[Tensor4 dtor] delete[] data: " << data_ << std::endl;
82 delete [] data_;
83 data_ = nullptr;
84 }
85 }
86
87 // Copy
88 Tensor4(const Tensor4& rhs)
89 {
90 if (rhs.ni_ <= 0 || rhs.nj_ <= 0 || rhs.nk_ <= 0 || rhs.nl_ <= 0) {
91 return;
92 }
93 allocate(rhs.ni_, rhs.nj_, rhs.nk_, rhs.nl_);
94 memcpy(data_, rhs.data_, size_*sizeof(T));
95 }
96
97 int num_i() { return ni_; }
98 int num_j() { return nj_; }
99 int num_k() { return nk_; }
100 int num_l() { return nl_; }
101
102 // Tensor4 assignment.
103 //
104 Tensor4& operator=(const Tensor4& rhs)
105 {
106 if (this == &rhs) {
107 return *this;
108 }
109
110 if (rhs.ni_ <= 0 || rhs.nj_ <= 0 || rhs.nk_ <= 0 || rhs.nl_ <= 0) {
111 return *this;
112 }
113
114 if (ni_ != rhs.ni_ || nj_ != rhs.nj_ || nk_ != rhs.nk_ || nl_ != rhs.nl_ <= 0) {
115 clear();
116 //delete [] data_;
117 //data_ = nullptr;
118 allocate(rhs.ni_, rhs.nj_, rhs.nk_, rhs.nl_);
119 }
120
121 memcpy(data_, rhs.data_, sizeof(T) * size_);
122 return *this;
123 }
124
125 Tensor4& operator=(const double value)
126 {
127 for (int i = 0; i < size_; i++) {
128 data_[i] = value;
129 }
130 return *this;
131 }
132
133 void write(const std::string& label, bool memory=true, T offset={}) const
134 {
135 int n = 0;
136 char file_prefix[1000];
137 for (auto c : label) {
138 if (c == '[') {
139 continue;
140 }
141 if (c == ']') {
142 file_prefix[n] = '_';
143 n += 1;
144 continue;
145 }
146 if (c == ' ') {
147 continue;
148 }
149 if (c == ':') {
150 c = '_';
151 }
152 file_prefix[n] = c;
153 n += 1;
154 }
155 file_prefix[n] = '\0';
156
157 FILE* fp;
158 std::string file_name;
159
160 // Write binary file.
161 file_name = std::string(file_prefix) + "_cm.bin";
162 std::cout << "[Tensor4::write] file_name: " << file_name << std::endl;
163 fp = fopen(file_name.c_str(), "wb");
164 fwrite(&size_, sizeof(int), 1, fp);
165 fwrite(data_, sizeof(double), size_, fp);
166 fclose(fp);
167 }
168
169 friend std::ostream& operator << (std::ostream& out, const Tensor4<T>& lhs)
170 {
171 for (int i = 0; i < lhs.size_; i++) {
172 out << lhs.data_[i];
173 if (i != lhs.size_-1) {
174 out << ", ";
175 }
176 }
177 return out;
178 }
179
180 //----------
181 // allocate
182 //----------
183 //
184 void allocate(const int num_i, const int num_j, const int num_k, const int num_l)
185 {
186 //std::cout << "----- Tensor4::allocate -----" << std::endl;
187 ni_ = num_i;
188 nj_ = num_j;
189 nk_ = num_k;
190 nl_ = num_l;
191 p1_ = num_i * num_j;
192 p2_ = p1_ * num_l;
193 size_ = ni_ * nj_ * nk_ * nl_;
194 data_ = new T [size_];
195 memset(data_, 0, sizeof(T)*size_);
196 //std::cout << "[Tensor4::allocate] data_: " << data_ << std::endl;
197 }
198
199 //-------------
200 // check_index
201 //-------------
202 //
203 void check_index(const int i, const int j, const int k, const int l) const
204 {
205 if (show_index_check_message) {
206 std::cout << "[Tensor4] **********************************" << std::endl;
207 std::cout << "[Tensor4] WARNING: Index checking is enabled " << std::endl;
208 std::cout << "[Tensor4] **********************************" << std::endl;
209 show_index_check_message = false;
210 }
211
212 if (data_ == nullptr) {
213 throw std::runtime_error(name_+"Accessing null data in Tensor4.");
214 }
215 if ((i < 0) || (i >= ni_) || (j < 0) || (j >= nj_) || (k < 0) || (k >= nk_) || (l < 0) || (l >= nl_)) {
216 auto i_str = std::to_string(ni_);
217 auto j_str = std::to_string(nj_);
218 auto k_str = std::to_string(nk_);
219 auto l_str = std::to_string(nl_);
220 auto dims = i_str + " x " + j_str + " x " + k_str + " x " + l_str;
221 auto index_str = " " + std::to_string(i) + "," + std::to_string(j) + "," + std::to_string(k) + "," + std::to_string(l);
222 throw std::runtime_error("Index (i,j,k,l)=" + index_str + " is out of bounds for " + dims + " array.");
223 }
224 }
225
226 //-------
227 // clear
228 //-------
229 // Free the array data.
230 //
231 // This is to replicate the Fortran DEALLOCATE().
232 //
233 void clear()
234 {
235 //std::cout << "----- Tensor4::erase -----" << std::endl;
236 if (data_ != nullptr) {
237 //std::cout << "[Tensor4::erase] data_: " << data_ << std::endl;
238 delete [] data_;
239 }
240
241 ni_ = 0;
242 nj_ = 0;
243 nk_ = 0;
244 nl_ = 0;
245 p1_ = 0;
246 p2_ = 0;
247 size_ = 0;
248 data_ = nullptr;
249 }
250
251 //--------
252 // resize
253 //--------
254 // Resize the array.
255 //
256 void resize(const int num_i, const int num_j, const int num_k, const int num_l)
257 {
258 if (data_ != nullptr) {
259 //std::cout << "[Tensor4::resize] data_: " << data_ << std::endl;
260 delete [] data_;
261 data_ = nullptr;
262 }
263
264 allocate(num_i, num_j, num_k, num_l);
265 }
266
267 /////////////////////////
268 // O p e r a t o r s //
269 /////////////////////////
270
271 const T& operator()(const int i, const int j, const int k, const int l) const
272 {
273 #ifdef Tensor4_check_enabled
274 check_index(i, j, k , l);
275 #endif
276 return data_[i + j*ni_ + p1_*k + p2_*l ];
277 }
278
279 T& operator()(const int i, const int j, const int k, const int l)
280 {
281 #ifdef Tensor4_check_enabled
282 check_index(i, j, k , l);
283 #endif
284 return data_[i + j*ni_ + p1_*k + p2_*l ];
285 }
286
287 ///////////////////////////////////////////////////
288 // M a t h e m a t i c a l o p e r a t o r s //
289 ///////////////////////////////////////////////////
290
291 // Multiply by a scalar.
292 //
293 Tensor4<T> operator*(const T value) const
294 {
295 Tensor4<T> result(ni_, nj_, nk_, nl_);
296 for (int i = 0; i < size_; i++) {
297 result.data_[i] = value * data_[i];
298 }
299 return result;
300 }
301
302 friend const Tensor4<T> operator * (const T value, const Tensor4& rhs)
303 {
304 if (rhs.data_ == nullptr) {
305 throw std::runtime_error("Null data for rhs Tensor4.");
306 }
307 Tensor4<T> result(rhs.ni_, rhs.nj_, rhs.nk_, rhs.nl_);
308 for (int i = 0; i < rhs.size_; i++) {
309 result.data_[i] = value * rhs.data_[i];
310 }
311 return result;
312 }
313
314 // Compound add assignment.
315 //
316 Tensor4<T> operator+=(const Tensor4<T>& rhs) const
317 {
318 for (int i = 0; i < size_; i++) {
319 data_[i] += rhs.data_[i];
320 }
321 return *this;
322 }
323
324 // Compound subtract assignment.
325 //
326 Tensor4<T> operator-=(const Tensor4<T>& rhs) const
327 {
328 for (int i = 0; i < size_; i++) {
329 data_[i] -= rhs.data_[i];
330 }
331 return *this;
332 }
333
334 // Add and subtract arrays.
335 //
336 Tensor4<T> operator+(const Tensor4<T>& rhs) const
337 {
338 Tensor4<T> result(rhs.ni_, rhs.nj_, rhs.nk_, rhs.nl_);
339 for (int i = 0; i < size_; i++) {
340 result.data_[i] = data_[i] + rhs.data_[i];
341 }
342 return result;
343 }
344
345 Tensor4<T> operator-(const Tensor4<T>& rhs) const
346 {
347 Tensor4<T> result(rhs.ni_, rhs.nj_, rhs.nk_, rhs.nl_);
348 for (int i = 0; i < size_; i++) {
349 result.data_[i] = data_[i] - rhs.data_[i];
350 }
351 return result;
352 }
353
354};
355
356#endif
357
The Tensor4 template class implements a simple interface to 4th order tensors.
Definition Tensor4.h:45