Multiscale Universal Interface  2.0
A Concurrent Framework for Coupling Heterogeneous Solvers
sampler_rbf.h
Go to the documentation of this file.
1 /*****************************************************************************
2  * Multiscale Universal Interface Code Coupling Library *
3  * *
4  * Copyright (C) 2019 Y. H. Tang, S. Kudo, X. Bian, Z. Li, G. E. Karniadakis, *
5  * A. Skillen, W. Liu, S. Longshaw, O. Mahfoze *
6  * *
7  * This software is jointly licensed under the Apache License, Version 2.0 *
8  * and the GNU General Public License version 3, you may use it according *
9  * to either. *
10  * *
11  * *
12  * Licensed under the Apache License, Version 2.0 (the "License"); *
13  * you may not use this file except in compliance with the License. *
14  * You may obtain a copy of the License at *
15  * *
16  * http://www.apache.org/licenses/LICENSE-2.0 *
17  * *
18  * Unless required by applicable law or agreed to in writing, software *
19  * distributed under the License is distributed on an "AS IS" BASIS, *
20  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
21  * See the License for the specific language governing permissions and *
22  * limitations under the License. *
23  * *
24  * ** GNU General Public License, version 3 ** *
25  * *
26  * This program is free software: you can redistribute it and/or modify *
27  * it under the terms of the GNU General Public License as published by *
28  * the Free Software Foundation, either version 3 of the License, or *
29  * (at your option) any later version. *
30  * *
31  * This program is distributed in the hope that it will be useful, *
32  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
33  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
34  * GNU General Public License for more details. *
35  * *
36  * You should have received a copy of the GNU General Public License *
37  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
38  *****************************************************************************/
39 
47 #ifndef MUI_SAMPLER_RBF_H_
48 #define MUI_SAMPLER_RBF_H_
49 
50 #include "../../general/util.h"
51 #include "../../uniface.h"
52 #include "../../linear_algebra/solver.h"
53 #include <iterator>
54 #include <ctime>
55 #include <sys/types.h>
56 #include <sys/stat.h>
57 #include <mpi.h>
58 
59 namespace mui {
60 
61 template<typename CONFIG = default_config,
62 typename O_TP = typename CONFIG::REAL, typename I_TP = O_TP>
63 class sampler_rbf {
64 
65 public:
66  using OTYPE = O_TP;
67  using ITYPE = I_TP;
68  using REAL = typename CONFIG::REAL;
69  using INT = typename CONFIG::INT;
70  using point_type = typename CONFIG::point_type;
71  using EXCEPTION = typename CONFIG::EXCEPTION;
72 
73  static const bool QUIET = CONFIG::QUIET;
74  static const bool DEBUG = CONFIG::DEBUG;
75 
126  sampler_rbf(REAL r, const std::vector<point_type> &pts, INT basisFunc = 0,
127  bool conservative = false, bool smoothFunc = false, bool generateMatrix = true,
128  const std::string &writeFileAddress = std::string(), REAL cutOff = 1e-9,
129  REAL cgSolveTol = 1e-6, INT cgMaxIter = 0, INT pouSize = 50,
130  INT precond = 1, MPI_Comm local_comm = MPI_COMM_NULL) :
131  r_(r),
132  pts_(pts),
133  basisFunc_(basisFunc),
134  conservative_(conservative),
135  consistent_(!conservative),
136  smoothFunc_(smoothFunc),
137  generateMatrix_(generateMatrix),
138  writeFileAddress_(writeFileAddress),
139  precond_(precond),
140  initialised_(false),
141  local_mpi_comm_world_(local_comm),
142  CABrow_(0),
143  CABcol_(0),
144  Hrow_(0),
145  Hcol_(0),
146  pouEnabled_(pouSize == 0 ? false : true),
147  cgSolveTol_(cgSolveTol),
148  cgMaxIter_(cgMaxIter),
149  N_sp_(pouSize),
150  M_ap_(pouSize),
151  local_rank_(0),
152  local_size_(0) {
153  //set s to give rbf(r)=cutOff (default 1e-9)
154  s_ = std::pow(-std::log(cutOff), 0.5) / r_;
155  twor_ = r_ * r_;
156  // Ensure CG solver parameters are sensible
157  if (cgMaxIter_ < 0)
158  cgMaxIter_ = 0;
159  if (cgSolveTol_ < 0)
160  cgSolveTol_ = 0;
161 
162  // Initialise the extened local points (ptsExtend_) by assign local points (pts_)
163  ptsExtend_.assign(pts_.begin(), pts_.end());
164 
165  if (local_mpi_comm_world_ != MPI_COMM_NULL) {
166  MPI_Comm_size(local_mpi_comm_world_, &local_size_);
167  MPI_Comm_rank(local_mpi_comm_world_, &local_rank_);
168  }
169 
170  // Create the matrix output folder if a directory name is given
171  if ( !writeFileAddress_.empty() ) {
172  int createError = 0;
173  #if defined(_WIN32) //Handle creation in Windows
174  createError = _mkdir(writeFileAddress_.c_str());
175  #else //Handle creation otherwise
176  createError = mkdir(writeFileAddress_.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
177  #endif
178  if (createError != 0 && createError != EEXIST)
179  EXCEPTION(std::runtime_error("MUI Error [sampler_rbf.h]: Problem creating RBF matrix folder"));
180  }
181  }
182 
183  template<template<typename, typename > class CONTAINER>
184  inline OTYPE filter(point_type focus, const CONTAINER<ITYPE, CONFIG> &data_points) const {
185  OTYPE sum = 0;
186 
187  // RBF matrix not yet created
188  if (!initialised_) {
189  if (generateMatrix_) { // Generating the matrix
190  const clock_t begin_time = clock();
191  facilitateGhostPoints();
192  REAL error = computeRBFtransformationMatrix(data_points, writeFileAddress_);
193 
194  if (!QUIET) {
195  std::cout << "MUI [sampler_rbf.h]: Matrices generated in: "
196  << static_cast<double>(clock() - begin_time) / CLOCKS_PER_SEC << "s ";
197  if (generateMatrix_) {
198  std::cout << std::endl
199  << " Average CG error: " << error << std::endl;
200  }
201  }
202  }
203  else // Matrix not found and not generating
204  EXCEPTION(std::runtime_error("MUI Error [sampler_rbf.h]: RBF matrix not found, call readRBFMatrix() first"));
205  }
206 
207  //Output for debugging
208  if ((!QUIET) && (DEBUG)) {
209  std::cout << "MUI [sampler_rbf.h]: Number of remote points: " << data_points.size()
210  << " at rank " << local_rank_ << " out of total ranks "
211  << local_size_ << std::endl;
212  std::cout << "MUI [sampler_rbf.h]: Number of local points: " << pts_.size()
213  << " at rank " << local_rank_ << " out of total ranks "
214  << local_size_ << std::endl;
215  std::cout << "MUI [sampler_rbf.h]: Number of ghost local points: " << ptsGhost_.size()
216  << " at rank " << local_rank_ << " out of total ranks "
217  << local_size_ << std::endl;
218  std::cout << "MUI [sampler_rbf.h]: Number of extended local points: "
219  << ptsExtend_.size() << " at rank " << local_rank_
220  << " out of total ranks " << local_size_ << std::endl;
221 
222  for (auto xPtsExtend : ptsExtend_) {
223  std::cout << " ";
224  int xPtsExtendSize;
225  try {
226  if (sizeof(xPtsExtend) == 0) {
227  throw "MUI Error [sampler_rbf.h]: Error zero xPtsExtend element exception";
228  }
229  else if (sizeof(xPtsExtend) < 0) {
230  throw "MUI Error [sampler_rbf.h]: Error invalid xPtsExtend element exception";
231  }
232  else if (sizeof(xPtsExtend[0]) == 0) {
233  throw "MUI Error [sampler_rbf.h]: Division by zero value of xPtsExtend[0] exception";
234  }
235  else if (sizeof(xPtsExtend[0]) < 0) {
236  throw "MUI Error [sampler_rbf.h]: Division by invalid value of xPtsExtend[0] exception";
237  }
238  xPtsExtendSize = sizeof(xPtsExtend) / sizeof(xPtsExtend[0]);
239  }
240  catch (const char *msg) {
241  std::cerr << msg << std::endl;
242  }
243  for (int i = 0; i < xPtsExtendSize; ++i) {
244  std::cout << xPtsExtend[i] << " ";
245  }
246  std::cout << std::endl;
247  }
248  }
249 
250  auto p = std::find_if(pts_.begin(), pts_.end(), [focus](point_type b) {
251  return normsq(focus - b) < std::numeric_limits<REAL>::epsilon();
252  });
253 
254  if (p == std::end(pts_))
255  EXCEPTION(std::runtime_error("MUI Error [sampler_rbf.h]: Point not found. Must pre-set points for RBF interpolation"));
256 
257  auto i = std::distance(pts_.begin(), p);
258  REAL tolerance = 1e-5;
259 
260  for (size_t j = 0; j < data_points.size(); j++) {
265  OTYPE HElement = 0;
266 
267  if (normsq(remote_pts_[j] - data_points[j].first) < (std::numeric_limits<REAL>::epsilon() + tolerance)){
268  HElement = H_.get_value(i, j);
269  }
270  else {
271 
272  point_type data_p = data_points[j].first;
273 
274  auto remote_p = std::find_if(remote_pts_.begin(), remote_pts_.end(), [data_p, tolerance](point_type b) {
275  return normsq(data_p - b) < (std::numeric_limits<REAL>::epsilon() + tolerance);
276  });
277 
278  if (remote_p == std::end(remote_pts_)) {
279  std::cout<< "Missing Remote points: " << data_points[j].first[0] << " " << data_points[j].first[1] << " Size of remote_pts_: " << remote_pts_.size()<< " at rank: " << local_rank_ << std::endl;
280  EXCEPTION(std::runtime_error("MUI Error [sampler_rbf.h]: Remote point not found. Must use the same set of remote points as used to construct the RBF coupling matrix"));
281  }
282 
283  auto new_j = std::distance(remote_pts_.begin(), remote_p);
284 
285  HElement = H_.get_value(i, new_j);
286  }
287 
288  sum += HElement * data_points[j].second;
289  }
290 
291  return sum;
292  }
293 
294  inline geometry::any_shape<CONFIG> support(point_type focus, REAL domain_mag) const {
295  return geometry::any_shape<CONFIG>();
296  }
297 
298  inline void preSetFetchPoints(std::vector<point_type> &pts) {
299  pts_ = pts;
300  initialised_ = false;
301  }
302 
303  inline void preSetFetchPointsExtend(std::vector<point_type> &pts) {
304  ptsExtend_ = pts;
305  initialised_ = false;
306  }
307 
308  inline void addFetchPoint(point_type pt) {
309  pts_.emplace_back(pt);
310  initialised_ = false;
311  }
312 
313  inline void addFetchPointExtend(point_type pt) {
314  ptsExtend_.emplace_back(pt);
315  initialised_ = false;
316  }
317 
318  inline void addFetchPointGhost(point_type pt) {
319  ptsGhost_.emplace_back(pt);
320  initialised_ = false;
321  }
322 
323  inline void readRBFMatrix(const std::string& readFileAddress) const {
324  std::ifstream inputFileMatrixSize(readFileAddress + "/matrixSize.dat");
325 
326  if (!inputFileMatrixSize) {
327  std::cerr << "MUI Error [sampler_rbf.h]: Could not locate the file address of matrixSize.dat"
328  << std::endl;
329  }
330  else {
331  std::string tempS;
332  std::vector<INT> tempV;
333  while (std::getline(inputFileMatrixSize, tempS)) {
334  // Skips the line if the first two characters are '//'
335  if (tempS[0] == '/' && tempS[1] == '/')
336  continue;
337  std::stringstream lineStream(tempS);
338  std::string tempSS;
339  while (std::getline(lineStream, tempSS, ',')) {
340  tempV.emplace_back(std::stoi(tempSS));
341  }
342  }
343  CABrow_ = tempV[0];
344  CABcol_ = tempV[1];
345  CAArow_ = tempV[2];
346  CAAcol_ = tempV[3];
347  Hrow_ = tempV[4];
348  Hcol_ = tempV[5];
349  remote_pts_num_ = tempV[6];
350  remote_pts_dim_ = tempV[7];
351  }
352 
353  std::ifstream inputFileCAB(readFileAddress + "/connectivityAB.dat");
354 
355  if (!inputFileCAB) {
356  std::cerr
357  << "MUI Error [sampler_rbf.h]: Could not locate the file address on the connectivityAB.dat"
358  << std::endl;
359  }
360  else {
361  connectivityAB_.resize(CABrow_);
362  for (INT i = 0; i < CABrow_; i++) {
363  connectivityAB_[i].resize(CABcol_, -1);
364  std::string tempS;
365  while (std::getline(inputFileCAB, tempS)) {
366  // Skips the line if the first two characters are '//'
367  if (tempS[0] == '/' && tempS[1] == '/')
368  continue;
369  std::stringstream lineStream(tempS);
370  std::string tempSS;
371  std::vector<INT> tempV;
372  while (std::getline(lineStream, tempSS, ',')) {
373  tempV.emplace_back(std::stoi(tempSS));
374  }
375  connectivityAB_.emplace_back(tempV);
376  }
377  }
378  }
379 
380  if (smoothFunc_) {
381  std::ifstream inputFileCAA(readFileAddress + "/connectivityAA.dat");
382 
383  if (!inputFileCAA) {
384  std::cerr
385  << "MUI Error [sampler_rbf.h]: Could not locate the file address on the connectivityAA.dat"
386  << std::endl;
387  }
388  else {
389  if ((CAArow_ == 0) || (CAAcol_ == 0)) {
390  std::cerr
391  << "MUI Error [sampler_rbf.h]: Error on the size of connectivityAA matrix in matrixSize.dat. Number of rows: "
392  << CAArow_ << " number of columns: " << CAAcol_
393  << ". Make sure matrices were generated with the smoothing function switched on."
394  << std::endl;
395  }
396  else {
397  connectivityAA_.resize(CAArow_);
398 
399  for (INT i = 0; i < CAArow_; i++) {
400  connectivityAA_[i].resize(CAAcol_, -1);
401  std::string tempS;
402  while (std::getline(inputFileCAA, tempS)) {
403  // Skips the line if the first two characters are '//'
404  if (tempS[0] == '/' && tempS[1] == '/')
405  continue;
406  std::stringstream lineStream(tempS);
407  std::string tempSS;
408  std::vector<INT> tempV;
409  while (std::getline(lineStream, tempSS, ',')) {
410  tempV.emplace_back(std::stoi(tempSS));
411  }
412  connectivityAA_.emplace_back(tempV);
413  }
414  }
415  }
416  }
417  }
418 
419  H_.resize(Hrow_, Hcol_);
420  H_.set_zero();
421 
422  std::ifstream inputFileHMatrix(readFileAddress + "/Hmatrix.dat");
423 
424  if (!inputFileHMatrix) {
425  std::cerr << "MUI Error [sampler_rbf.h]: Could not locate the file address on the Hmatrix.dat"
426  << std::endl;
427  }
428  else {
429  inputFileHMatrix >> H_;
430 
431  if ((H_.get_rows() != Hrow_) || (H_.get_cols() != Hcol_)) {
432  std::cerr << "row of H_ (" << H_.get_rows()
433  << ") is not NOT equal to Hrow_ (" << Hrow_ << "), or"
434  << std::endl << "column of H_ (" << H_.get_cols()
435  << ") is not NOT equal to Hcol_ (" << Hcol_ << ")"
436  << std::endl;
437  }
438  }
439 
440  std::ifstream inputFileRemotePoints(readFileAddress + "/remotePoints.dat");
441 
442  if (!inputFileRemotePoints) {
443  std::cerr << "MUI Error [sampler_rbf.h]: Could not locate the file address on the Hmatrix.dat"
444  << std::endl;
445  }
446  else {
447 
448  if (CONFIG::D != remote_pts_dim_){
449  EXCEPTION(std::runtime_error("MUI Error [sampler_rbf.h]: CONFIG::D must equal to remote point dimension in remotePoints.dat"));
450  }
451 
452  std::string tempS;
453  while (std::getline(inputFileRemotePoints, tempS)) {
454  // Skips the line if the first two characters are '//'
455  if (tempS[0] == '/' && tempS[1] == '/')
456  continue;
457  std::stringstream lineStream(tempS);
458  std::string tempSS;
459  point_type tempP;
460  size_t temp_i = 0;
461  while (std::getline(lineStream, tempSS, ',')) {
462  assert(temp_i<static_cast<size_t>(remote_pts_dim_));
463  tempP[temp_i] = static_cast<REAL>(std::stold(tempSS));
464  temp_i++;
465  }
466  remote_pts_.emplace_back(tempP);
467  }
468 
469  assert(remote_pts_.size() == static_cast<size_t>(remote_pts_num_));
470  }
471 
472  initialised_ = true;
473  }
474 
475  // Functions to facilitate ghost points
476 
477  // Determine bounding box of local points
478  std::pair<point_type, point_type> localBoundingBox(const std::vector<point_type> pt) const {
479  point_type lbbMin, lbbMax;
480  REAL maxVal = std::numeric_limits<REAL>::max();
481  try {
482  if (CONFIG::D == 1) {
483  lbbMax[0] = -maxVal;
484  lbbMin[0] = maxVal;
485  }
486  else if (CONFIG::D == 2) {
487  lbbMax[0] = -maxVal;
488  lbbMax[1] = -maxVal;
489  lbbMin[0] = maxVal;
490  lbbMin[1] = maxVal;
491  }
492  else if (CONFIG::D == 3) {
493  lbbMax[0] = -maxVal;
494  lbbMax[1] = -maxVal;
495  lbbMax[2] = -maxVal;
496  lbbMin[0] = maxVal;
497  lbbMin[1] = maxVal;
498  lbbMin[2] = maxVal;
499  }
500  else {
501  throw "MUI Error [sampler_rbf.h]: Invalid value of CONFIG::D exception";
502  }
503  }
504  catch (const char *msg) {
505  std::cerr << msg << std::endl;
506  }
507 
508  for (auto xPts : pt) {
509  try {
510  switch (CONFIG::D) {
511  case 1:
512  if (xPts[0] > lbbMax[0])
513  lbbMax[0] = xPts[0];
514  if (xPts[0] < lbbMin[0])
515  lbbMin[0] = xPts[0];
516  break;
517  case 2:
518  if (xPts[0] > lbbMax[0])
519  lbbMax[0] = xPts[0];
520  if (xPts[0] < lbbMin[0])
521  lbbMin[0] = xPts[0];
522  if (xPts[1] > lbbMax[1])
523  lbbMax[1] = xPts[1];
524  if (xPts[1] < lbbMin[1])
525  lbbMin[1] = xPts[1];
526  break;
527  case 3:
528  if (xPts[0] > lbbMax[0])
529  lbbMax[0] = xPts[0];
530  if (xPts[0] < lbbMin[0])
531  lbbMin[0] = xPts[0];
532  if (xPts[1] > lbbMax[1])
533  lbbMax[1] = xPts[1];
534  if (xPts[1] < lbbMin[1])
535  lbbMin[1] = xPts[1];
536  if (xPts[2] > lbbMax[2])
537  lbbMax[2] = xPts[2];
538  if (xPts[2] < lbbMin[2])
539  lbbMin[2] = xPts[2];
540  break;
541  default:
542  throw "MUI Error [sampler_rbf.h]: Invalid value of CONFIG::D exception";
543  }
544  }
545  catch (const char *msg) {
546  std::cerr << msg << std::endl;
547  }
548  }
549  return std::make_pair(lbbMin, lbbMax);
550  }
551 
552  // Determine extended bounding box of local points include ghost area
553  std::pair<point_type, point_type> localExtendBoundingBox(std::pair<point_type, point_type> lbb, REAL r) const {
554  point_type lbbExtendMin, lbbExtendMax;
555  lbbExtendMin = lbb.first;
556  lbbExtendMax = lbb.second;
557  for (INT i = 0; i < CONFIG::D; ++i) {
558  lbbExtendMax[i] += r;
559  lbbExtendMin[i] -= r;
560  }
561  return std::make_pair(lbbExtendMin, lbbExtendMax);
562  }
563 
564  // Collect points that will be sent to other processors as ghost points
565  std::pair<std::vector<std::pair<INT, INT>>, std::vector<std::pair<INT, std::vector<point_type>>>> getGhostPointsToSend(
566  std::pair<point_type, point_type> lbbExtend, MPI_Comm local_world,
567  int local_rank, int local_size) const {
568  std::pair<INT, std::pair<point_type, point_type>> localBB, globalBB[local_size];
569 
570  // Assign rank ID and extended bounding box of local points to localBB
571  localBB.first = local_rank;
572  localBB.second.first = lbbExtend.second;
573  localBB.second.second = lbbExtend.first;
574 
575  // Gather the bounding boxes from all processes
576  MPI_Allgather(&localBB,
577  sizeof(std::pair<INT, std::pair<point_type, point_type>>),
578  MPI_BYTE, globalBB,
579  sizeof(std::pair<INT, std::pair<point_type, point_type>>),
580  MPI_BYTE, local_world);
581 
582  // output for debugging
583  if ((!QUIET) && (DEBUG)) {
584  for (auto xGlobalBB : globalBB) {
585  std::cout << "MUI [sampler_rbf.h]: globalBBs: " << xGlobalBB.first << " "
586  << xGlobalBB.second.first[0] << " "
587  << xGlobalBB.second.first[1] << " "
588  << xGlobalBB.second.second[0] << " "
589  << xGlobalBB.second.second[1] << " at rank "
590  << local_rank_ << std::endl;
591  }
592  }
593 
594  // Declare vector to gather number of points to send to each processors
595  std::vector<std::pair<INT, INT>> ghostPointsCountToSend;
596  for (auto xGlobalBB : globalBB) {
597  if (xGlobalBB.first == local_rank)
598  continue;
599  ghostPointsCountToSend.push_back(std::make_pair(xGlobalBB.first, 0));
600  }
601 
602  // Loop over local points to collect ghost points for other processors
603  std::vector<std::pair<INT, std::vector<point_type>>> ghostPointsToSend;
604  for (size_t i = 0; i < pts_.size(); ++i) {
605  for (auto xGlobalBB : globalBB) {
606  if (xGlobalBB.first == local_rank)
607  continue;
608  try {
609  switch (CONFIG::D) {
610  case 1: {
611  if ((pts_[i][0] <= xGlobalBB.second.first[0])
612  && (pts_[i][0] >= xGlobalBB.second.second[0])) {
613 
614  auto ghostPointsToSendIter =
615  std::find_if(ghostPointsToSend.begin(),
616  ghostPointsToSend.end(),
617  [xGlobalBB](std::pair<INT, std::vector<point_type>> b) {
618  return b.first == xGlobalBB.first;
619  });
620 
621  auto ghostPointsCountToSendIter = std::find_if(
622  ghostPointsCountToSend.begin(),
623  ghostPointsCountToSend.end(),
624  [xGlobalBB](std::pair<INT, INT> b) {
625  return b.first == xGlobalBB.first;
626  });
627 
628  if (ghostPointsToSendIter
629  == std::end(ghostPointsToSend)) {
630  std::vector<point_type> vecPointTemp { pts_[i] };
631  ghostPointsToSend.push_back(std::make_pair(xGlobalBB.first, vecPointTemp));
632  }
633  else
634  ghostPointsToSendIter->second.push_back(pts_[i]);
635 
636  assert(ghostPointsCountToSendIter != std::end(ghostPointsCountToSend));
637 
638  ++ghostPointsCountToSendIter->second;
639  }
640  break;
641  }
642  case 2: {
643  if ((pts_[i][0] <= xGlobalBB.second.first[0])
644  && (pts_[i][0] >= xGlobalBB.second.second[0])
645  && (pts_[i][1] <= xGlobalBB.second.first[1])
646  && (pts_[i][1] >= xGlobalBB.second.second[1])) {
647 
648  auto ghostPointsToSendIter =
649  std::find_if(ghostPointsToSend.begin(),
650  ghostPointsToSend.end(),
651  [xGlobalBB](std::pair<INT, std::vector<point_type>> b) {
652  return b.first == xGlobalBB.first;
653  });
654 
655  auto ghostPointsCountToSendIter = std::find_if(
656  ghostPointsCountToSend.begin(),
657  ghostPointsCountToSend.end(),
658  [xGlobalBB](std::pair<INT, INT> b) {
659  return b.first == xGlobalBB.first;
660  });
661 
662  if (ghostPointsToSendIter
663  == std::end(ghostPointsToSend)) {
664  std::vector<point_type> vecPointTemp { pts_[i] };
665  ghostPointsToSend.push_back(std::make_pair(xGlobalBB.first, vecPointTemp));
666  }
667  else
668  ghostPointsToSendIter->second.push_back(pts_[i]);
669 
670  assert(ghostPointsCountToSendIter != std::end(ghostPointsCountToSend));
671 
672  ++ghostPointsCountToSendIter->second;
673  }
674  break;
675  }
676  case 3: {
677  if ((pts_[i][0] <= xGlobalBB.second.first[0])
678  && (pts_[i][0] >= xGlobalBB.second.second[0])
679  && (pts_[i][1] <= xGlobalBB.second.first[1])
680  && (pts_[i][1] >= xGlobalBB.second.second[1])
681  && (pts_[i][2] <= xGlobalBB.second.first[2])
682  && (pts_[i][2] >= xGlobalBB.second.second[2])) {
683 
684  auto ghostPointsToSendIter =
685  std::find_if(ghostPointsToSend.begin(),
686  ghostPointsToSend.end(),
687  [xGlobalBB](std::pair<INT, std::vector<point_type>> b) {
688  return b.first == xGlobalBB.first;
689  });
690 
691  auto ghostPointsCountToSendIter = std::find_if(
692  ghostPointsCountToSend.begin(),
693  ghostPointsCountToSend.end(),
694  [xGlobalBB](std::pair<INT, INT> b) {
695  return b.first == xGlobalBB.first;
696  });
697 
698  if (ghostPointsToSendIter
699  == std::end(ghostPointsToSend)) {
700  std::vector<point_type> vecPointTemp { pts_[i] };
701  ghostPointsToSend.push_back(std::make_pair(xGlobalBB.first, vecPointTemp));
702  }
703  else
704  ghostPointsToSendIter->second.push_back(pts_[i]);
705 
706  assert(ghostPointsCountToSendIter != std::end(ghostPointsCountToSend));
707 
708  ++ghostPointsCountToSendIter->second;
709  }
710  break;
711  }
712  default:
713  throw "MUI Error [sampler_rbf.h]: Invalid value of CONFIG::D exception";
714  }
715  }
716  catch (const char *msg) {
717  std::cerr << msg << std::endl;
718  }
719  }
720  }
721 
722  // output for debugging
723  if ((!QUIET) && (DEBUG)) {
724  std::cout << "MUI [sampler_rbf.h]: Total size of GhostPointsToSend "
725  << ghostPointsToSend.size() << " at rank " << local_rank
726  << std::endl;
727  for (auto xGhostPointsToSend : ghostPointsToSend) {
728  std::cout << "MUI [sampler_rbf.h]: xGhostPointsToSend to rank "
729  << xGhostPointsToSend.first << " at rank " << local_rank
730  << std::endl;
731  for (auto xVectPts : xGhostPointsToSend.second) {
732  std::cout << "MUI [sampler_rbf.h]: " << xVectPts[0] << " " << xVectPts[1]
733  << " at rank " << local_rank << std::endl;
734  }
735  }
736  for (auto xGhostPointsCountToSend : ghostPointsCountToSend) {
737  std::cout << "MUI [sampler_rbf.h]: xGhostPointsCountToSend to rank "
738  << xGhostPointsCountToSend.first << " has "
739  << xGhostPointsCountToSend.second
740  << " points to send at rank " << local_rank_
741  << std::endl;
742  }
743  }
744  return std::make_pair(ghostPointsCountToSend, ghostPointsToSend);
745  }
746 
747  // Distribution of ghost points among processors by all to all
748  std::vector<point_type> distributeGhostPoints(
749  std::vector<std::pair<INT, INT>> ghostPointsCountToSend,
750  std::vector<std::pair<INT, std::vector<point_type>>> ghostPointsToSend,
751  MPI_Comm local_world, int local_rank) const {
752  std::vector<point_type> ptsGhost;
753  std::vector<std::pair<INT, INT>> ghostPointsCountToRecv;
754  for (auto xGhostPointsCountToSend : ghostPointsCountToSend) {
755  assert(xGhostPointsCountToSend.first != local_rank);
756  assert(xGhostPointsCountToSend.second >= 0);
757 
758  // Determined of number of points to transfer by pairwise communication
759  MPI_Send(&xGhostPointsCountToSend.second, 1, MPI_INT,
760  xGhostPointsCountToSend.first, 0, local_world);
761  int pointsCountTemp = -1;
762  MPI_Recv(&pointsCountTemp, 1, MPI_INT,
763  xGhostPointsCountToSend.first, 0, local_world,
764  MPI_STATUS_IGNORE);
765  assert(pointsCountTemp >= 0);
766  ghostPointsCountToRecv.push_back(
767  std::make_pair(xGhostPointsCountToSend.first,
768  pointsCountTemp));
769 
770  // Send ghost points by pairwise communication
771  if (xGhostPointsCountToSend.second != 0) {
772  auto ghostPointsToSendIter = std::find_if(
773  ghostPointsToSend.begin(), ghostPointsToSend.end(),
774  [xGhostPointsCountToSend](
775  std::pair<INT, std::vector<point_type>> b) {
776  return b.first == xGhostPointsCountToSend.first;
777  });
778 
779  assert(ghostPointsToSendIter->second.size() == static_cast<size_t>(xGhostPointsCountToSend.second));
780 
781  int buffer_size = ghostPointsToSendIter->second.size() * sizeof(point_type);
782  char *buffer = new char[buffer_size];
783  int position = 0;
784  for (auto &pointElement : ghostPointsToSendIter->second) {
785  MPI_Pack(&pointElement, sizeof(point_type), MPI_BYTE,
786  buffer, buffer_size, &position, local_world);
787  }
788  MPI_Send(buffer, position, MPI_PACKED, xGhostPointsCountToSend.first, 1, local_world);
789  delete[] buffer;
790 
791  // output for debugging
792  if ((!QUIET) && (DEBUG)) {
793  for (auto xsend : ghostPointsToSendIter->second) {
794  std::cout << "MUI [sampler_rbf.h]: Send ghost point to rank "
795  << xGhostPointsCountToSend.first
796  << " with value of " << xsend[0] << " "
797  << xsend[1] << " at rank " << local_rank
798  << std::endl;
799  }
800  }
801  }
802 
803  // Receive ghost points by pairwise communication
804  auto ghostPointsCountToRecvIter = std::find_if(
805  ghostPointsCountToRecv.begin(),
806  ghostPointsCountToRecv.end(),
807  [xGhostPointsCountToSend](std::pair<INT, INT> b) {
808  return b.first == xGhostPointsCountToSend.first;
809  });
810 
811  if (ghostPointsCountToRecvIter->second != 0) {
812  std::vector<point_type> vecPointTemp;
813  MPI_Status status;
814  MPI_Probe(ghostPointsCountToRecvIter->first, 1, local_world,
815  &status);
816  int buffer_size;
817  MPI_Get_count(&status, MPI_PACKED, &buffer_size);
818  char *buffer = new char[buffer_size];
819  MPI_Recv(buffer, buffer_size, MPI_PACKED,
820  ghostPointsCountToRecvIter->first, 1, local_world,
821  MPI_STATUS_IGNORE);
822  int position = 0;
823  int count;
824  MPI_Get_elements(&status, MPI_BYTE, &count);
825  vecPointTemp.resize(count / sizeof(point_type));
826  for (auto &pointElement : vecPointTemp) {
827  MPI_Unpack(buffer, buffer_size, &position, &pointElement,
828  sizeof(point_type), MPI_BYTE, local_world);
829  }
830  delete[] buffer;
831  for (auto xvecPointTemp : vecPointTemp) {
832  ptsGhost.emplace_back(xvecPointTemp);
833  // output for debugging
834  if ((!QUIET) && (DEBUG)) {
835  std::cout << "MUI [sampler_rbf.h]: Receive ghost point from rank "
836  << ghostPointsCountToRecvIter->first
837  << " with value of " << xvecPointTemp[0] << " "
838  << xvecPointTemp[1] << " at rank " << local_rank
839  << std::endl;
840  }
841  }
842  }
843  }
844  return ptsGhost;
845  }
846 
847  // Facilitate Ghost points
848  void facilitateGhostPoints() const {
849  std::pair<point_type, point_type> lbb = localBoundingBox(pts_);
850  std::pair<point_type, point_type> lbbExtend = localExtendBoundingBox(lbb, r_);
851  if (local_mpi_comm_world_ != MPI_COMM_NULL) {
852  std::pair<std::vector<std::pair<INT, INT>>, std::vector<std::pair<INT, std::vector<point_type>>>>
853  ghostPointsToSendPair = getGhostPointsToSend(lbbExtend, local_mpi_comm_world_, local_rank_, local_size_);
854  ptsGhost_ = distributeGhostPoints(ghostPointsToSendPair.first,
855  ghostPointsToSendPair.second, local_mpi_comm_world_, local_rank_);
856  // output for debugging
857  if ((!QUIET) && (DEBUG)) {
858  std::cout << "MUI [sampler_rbf.h]: Local bounding box: " << lbb.second[0] << " "
859  << lbb.second[1] << " " << lbb.first[0] << " "
860  << lbb.first[1] << " at rank " << local_rank_
861  << " out of total ranks " << local_size_ << std::endl;
862  std::cout << "MUI [sampler_rbf.h]: Extended local bounding box: "
863  << lbbExtend.second[0] << " " << lbbExtend.second[1]
864  << " " << lbbExtend.first[0] << " "
865  << lbbExtend.first[1] << " at rank " << local_rank_
866  << " out of total ranks " << local_size_ << std::endl;
867  }
868  }
869  // Construct the extended local points by combine local points with ghost points
870  ptsExtend_.insert(ptsExtend_.end(), ptsGhost_.begin(), ptsGhost_.end());
871  }
872 
873 private:
874  template<template<typename, typename > class CONTAINER>
875  REAL computeRBFtransformationMatrix(const CONTAINER<ITYPE, CONFIG> &data_points, const std::string &fileAddress) const {
876  bool writeMatrix = !fileAddress.empty();
877 
878  // Refine partition size depending on if PoU enabled, whether conservative or consistent
879  // and if problem size smaller than defined patch size
880  if (conservative_) { // Conservative RBF, using local point set for size
881  if (pouEnabled_) { // PoU enabled so check patch size not larger than point set
882  if (ptsExtend_.size() < N_sp_)
883  N_sp_ = ptsExtend_.size();
884  }
885  else
886  N_sp_ = ptsExtend_.size();
887  }
888 
889  if (consistent_) { // Consistent RBF, using remote point set for size
890  if (pouEnabled_) { // PoU enabled so check patch size not larger than point set
891  if (data_points.size() < N_sp_)
892  N_sp_ = data_points.size();
893  }
894  else
895  N_sp_ = data_points.size();
896  }
897 
898  if (smoothFunc_) {
899  if (ptsExtend_.size() < M_ap_)
900  M_ap_ = ptsExtend_.size() - 1;
901  }
902 
903  REAL errorReturn = 0;
904 
905  if (conservative_)
906  buildConnectivityConservative(data_points, N_sp_, writeMatrix, fileAddress);
907  else
908  buildConnectivityConsistent(data_points, N_sp_, writeMatrix, fileAddress);
909 
910  H_.resize(ptsExtend_.size(), data_points.size());
911  H_.set_zero();
912 
913  if (smoothFunc_) {
914  buildConnectivitySmooth(M_ap_, writeMatrix, fileAddress);
915  H_toSmooth_.resize(ptsExtend_.size(), data_points.size());
916  H_toSmooth_.set_zero();
917  }
918 
919  if (writeMatrix) {
920  std::ofstream outputFileMatrixSize(fileAddress + "/matrixSize.dat");
921 
922  if (!outputFileMatrixSize) {
923  std::cerr
924  << "MUI Error [sampler_rbf.h]: Could not locate the file address of matrixSize.dat!"
925  << std::endl;
926  }
927  else {
928  outputFileMatrixSize
929  << "// *********************************************************************************************************************************************";
930  outputFileMatrixSize << "\n";
931  outputFileMatrixSize
932  << "// **** This is the 'matrixSize.dat' file of the RBF spatial sampler of the MUI library";
933  outputFileMatrixSize << "\n";
934  outputFileMatrixSize
935  << "// **** This file contains the size (number of rows and number of columns) of the Point Connectivity Matrix (N) and the Coupling Matrix (H).";
936  outputFileMatrixSize << "\n";
937  outputFileMatrixSize
938  << "// **** The file uses the Comma-Separated Values (CSV) format and the ASCII format with the meanings as follows: ";
939  outputFileMatrixSize << "\n";
940  outputFileMatrixSize
941  << "// **** The number of rows of the Point Connectivity Matrix (N), ";
942  outputFileMatrixSize << "\n";
943  outputFileMatrixSize
944  << "// **** The number of columns of the Point Connectivity Matrix (N),";
945  outputFileMatrixSize << "\n";
946  outputFileMatrixSize
947  << "// **** The number of rows of the Point Connectivity Matrix (M) (for smoothing), ";
948  outputFileMatrixSize << "\n";
949  outputFileMatrixSize
950  << "// **** The number of columns of the Point Connectivity Matrix (M) (for smoothing),";
951  outputFileMatrixSize << "\n";
952  outputFileMatrixSize
953  << "// **** The number of rows of the Coupling Matrix (H),";
954  outputFileMatrixSize << "\n";
955  outputFileMatrixSize
956  << "// **** The number of columns of the Coupling Matrix (H)";
957  outputFileMatrixSize << "\n";
958  outputFileMatrixSize
959  << "// **** The size of remote points";
960  outputFileMatrixSize << "\n";
961  outputFileMatrixSize
962  << "// **** The dimension of remote points";
963  outputFileMatrixSize << "\n";
964  outputFileMatrixSize
965  << "// *********************************************************************************************************************************************";
966  outputFileMatrixSize << "\n";
967  outputFileMatrixSize << "// ";
968  outputFileMatrixSize << "\n";
969  outputFileMatrixSize << connectivityAB_.size();
970  outputFileMatrixSize << ",";
971  outputFileMatrixSize << connectivityAB_[0].size();
972  outputFileMatrixSize << ",";
973  if (smoothFunc_) {
974  outputFileMatrixSize << connectivityAA_.size();
975  outputFileMatrixSize << ",";
976  outputFileMatrixSize << connectivityAA_[0].size();
977  outputFileMatrixSize << ",";
978  }
979  else {
980  outputFileMatrixSize << "0";
981  outputFileMatrixSize << ",";
982  outputFileMatrixSize << "0";
983  outputFileMatrixSize << ",";
984  }
985  outputFileMatrixSize << H_.get_rows();
986  outputFileMatrixSize << ",";
987  outputFileMatrixSize << H_.get_cols();
988  outputFileMatrixSize << ",";
989  outputFileMatrixSize << data_points.size();
990  outputFileMatrixSize << ",";
991  outputFileMatrixSize << CONFIG::D;
992  outputFileMatrixSize << "\n";
993  }
994  }
995 
996  if (conservative_) { // Build matrix for conservative RBF
997  errorReturn = buildMatrixConservative(data_points, N_sp_, M_ap_, smoothFunc_, pouEnabled_);
998  }
999  else {
1000  // Build matrix for consistent RBF
1001  errorReturn = buildMatrixConsistent(data_points, N_sp_, M_ap_, smoothFunc_, pouEnabled_);
1002  }
1003 
1004  if (writeMatrix) {
1005  std::ofstream outputFileHMatrix(fileAddress + "/Hmatrix.dat");
1006 
1007  if (!outputFileHMatrix) {
1008  std::cerr
1009  << "MUI Error [sampler_rbf.h]: Could not locate the file address of Hmatrix.dat!"
1010  << std::endl;
1011  }
1012  else {
1013  outputFileHMatrix
1014  << "// ************************************************************************************************";
1015  outputFileHMatrix << "\n";
1016  outputFileHMatrix
1017  << "// **** This is the 'Hmatrix.dat' file of the RBF spatial sampler of the MUI library";
1018  outputFileHMatrix << "\n";
1019  outputFileHMatrix
1020  << "// **** This file contains the entire matrix of the Coupling Matrix (H).";
1021  outputFileHMatrix << "\n";
1022  outputFileHMatrix
1023  << "// **** The file uses the Comma-Separated Values (CSV) format with ASCII for the entire H matrix";
1024  outputFileHMatrix << "\n";
1025  outputFileHMatrix
1026  << "// ************************************************************************************************";
1027  outputFileHMatrix << "\n";
1028  outputFileHMatrix << "// ";
1029  outputFileHMatrix << "\n";
1030  outputFileHMatrix << H_;
1031  }
1032  }
1033 
1034  initialised_ = true;
1035 
1036  return errorReturn;
1037  }
1038 
1039  template<template<typename, typename > class CONTAINER>
1040  inline REAL buildMatrixConsistent(const CONTAINER<ITYPE, CONFIG> &data_points, const size_t NP,
1041  const size_t MP, bool smoothing, bool pou) const {
1042  REAL errorReturn = 0;
1043  std::pair<INT, REAL> iterErrorReturn(0, 0);
1044  if( pou ) { // Using PoU approach
1045  for (size_t row = 0; row < ptsExtend_.size(); row++) {
1046  linalg::sparse_matrix<INT, REAL> Css; //< Matrix of radial basis function evaluations between prescribed points
1047  linalg::sparse_matrix<INT, REAL> Aas; //< Matrix of RBF evaluations between prescribed and interpolation points
1048 
1049  Css.resize((1 + NP + CONFIG::D), (1 + NP + CONFIG::D));
1050  Aas.resize((1 + NP + CONFIG::D), 1);
1051 
1052  // Set matrix Css
1053  // Define intermediate matrix for performance purpose
1054  linalg::sparse_matrix<INT, REAL> Css_coo((1 + NP + CONFIG::D),(1 + NP + CONFIG::D),"COO");
1055 
1056  for (size_t i = 0; i < NP; i++) {
1057  for (size_t j = i; j < NP; j++) {
1058  int glob_i = connectivityAB_[row][i];
1059  int glob_j = connectivityAB_[row][j];
1060 
1061  auto d = norm(
1062  data_points[glob_i].first
1063  - data_points[glob_j].first);
1064 
1065  if (d < r_) {
1066  REAL w = rbf(d);
1067  Css_coo.set_value(i, j, w, false);
1068 // Css.set_value(i, j, w);
1069  if (i != j) {
1070  Css_coo.set_value(j, i, w, false);
1071 // Css.set_value(j, i, w);
1072  }
1073  }
1074  }
1075  }
1076 
1077  for (size_t i = 0; i < NP; i++) {
1078  Css_coo.set_value(i, NP, 1, false);
1079  Css_coo.set_value(NP, i, 1, false);
1080 // Css.set_value(i, NP, 1);
1081 // Css.set_value(NP, i, 1);
1082 
1083  int glob_i = connectivityAB_[row][i];
1084 
1085  for (INT dim = 0; dim < CONFIG::D; dim++) {
1086  Css_coo.set_value(i, (NP + dim + 1),
1087  data_points[glob_i].first[dim], false);
1088 // Css.set_value(i, (NP + dim + 1),
1089 // data_points[glob_i].first[dim]);
1090  Css_coo.set_value((NP + dim + 1), i,
1091  data_points[glob_i].first[dim], false);
1092 // Css.set_value((NP + dim + 1), i,
1093 // data_points[glob_i].first[dim]);
1094  }
1095  }
1096  Css = Css_coo;
1097 
1098  // Set Aas
1099  // Define intermediate matrix for performance purpose
1100  linalg::sparse_matrix<INT, REAL> Aas_coo((1 + NP + CONFIG::D),1,"COO");
1101 
1102  for (size_t j = 0; j < NP; j++) {
1103  int glob_j = connectivityAB_[row][j];
1104 
1105  auto d = norm(ptsExtend_[row] - data_points[glob_j].first);
1106 
1107  if (d < r_) {
1108  Aas_coo.set_value(j, 0, rbf(d), false);
1109 // Aas.set_value(j, 0, rbf(d));
1110  }
1111  }
1112  Aas_coo.set_value(NP, 0, 1, false);
1113 // Aas.set_value(NP, 0, 1);
1114 
1115  for (int dim = 0; dim < CONFIG::D; dim++) {
1116  Aas_coo.set_value((NP + dim + 1), 0, ptsExtend_[row][dim], false);
1117 // Aas.set_value((NP + dim + 1), 0, ptsExtend_[row][dim]);
1118  }
1119  Aas = Aas_coo;
1120 
1121  linalg::sparse_matrix<INT, REAL> H_i;
1122 
1123  if (precond_ == 0) {
1124  linalg::conjugate_gradient<INT, REAL> cg(Css, Aas, cgSolveTol_, cgMaxIter_);
1125  iterErrorReturn = cg.solve();
1126  H_i = cg.getSolution();
1127  }
1128  else if (precond_ == 1) {
1129  linalg::diagonal_preconditioner<INT, REAL> M(Css);
1130  linalg::conjugate_gradient<INT, REAL> cg(Css, Aas, cgSolveTol_, cgMaxIter_, &M);
1131  iterErrorReturn = cg.solve();
1132  H_i = cg.getSolution();
1133  }
1134  else {
1135  std::cerr
1136  << "MUI Error [sampler_rbf.h]: Invalid CG Preconditioner function number ("
1137  << precond_ << ")" << std::endl
1138  << "Please set the CG Preconditioner function number (precond_) as: "
1139  << std::endl << "precond_=0 (No Preconditioner); "
1140  << std::endl << "precond_=1 (Diagonal Preconditioner); " << std::endl;
1141  std::abort();
1142  }
1143 
1144  if (DEBUG) {
1145  std::cout << "MUI [sampler_rbf.h]: #iterations of H_i: " << iterErrorReturn.first
1146  << ". Error of H_i: " << iterErrorReturn.second
1147  << std::endl;
1148  }
1149 
1150  errorReturn += iterErrorReturn.second;
1151 
1152  if (smoothing) {
1153  for (size_t j = 0; j < NP; j++) {
1154  INT glob_j = connectivityAB_[row][j];
1155  H_toSmooth_.set_value(row, glob_j, H_i.get_value(j, 0));
1156  }
1157  }
1158  else {
1159  for (size_t j = 0; j < NP; j++) {
1160  INT glob_j = connectivityAB_[row][j];
1161  H_.set_value(row, glob_j, H_i.get_value(j, 0));
1162  }
1163  }
1164  }
1165 
1166  errorReturn /= static_cast<REAL>(pts_.size());
1167 
1168  if (smoothing) {
1169  for (size_t row = 0; row < ptsExtend_.size(); row++) {
1170  for (size_t j = 0; j < NP; j++) {
1171  INT glob_j = connectivityAB_[row][j];
1172  REAL h_j_sum = 0.;
1173  REAL f_sum = 0.;
1174 
1175  for (size_t k = 0; k < MP; k++) {
1176  INT row_k = connectivityAA_[row][k];
1177  if (row_k == static_cast<INT>(row)) {
1178  std::cerr << "MUI Error [sampler_rbf.h]: Invalid row_k value: " << row_k
1179  << std::endl;
1180  }
1181  else
1182  h_j_sum += std::pow(dist_h_i(row, row_k), -2.);
1183  }
1184 
1185  for (size_t k = 0; k < MP; k++) {
1186  INT row_k = connectivityAA_[row][k];
1187  if (row_k == static_cast<INT>(row)) {
1188  std::cerr << "MUI Error [sampler_rbf.h]: Invalid row_k value: " << row_k
1189  << std::endl;
1190  }
1191  else {
1192  REAL w_i = ((std::pow(dist_h_i(row, row_k), -2.)) / (h_j_sum));
1193  f_sum += w_i * H_toSmooth_.get_value(row_k, glob_j);
1194  }
1195  }
1196  H_.set_value(row, glob_j, (0.5 * (f_sum + H_toSmooth_.get_value(row, glob_j))));
1197  }
1198  }
1199  }
1200  }
1201  else { // Not using PoU
1202  linalg::sparse_matrix<INT,REAL> Css; //< Matrix of radial basis function evaluations between prescribed points
1203  linalg::sparse_matrix<INT,REAL> Aas; //< Matrix of RBF evaluations between prescribed and interpolation points
1204 
1205  Css.resize((1 + data_points.size() + CONFIG::D), (1 + data_points.size() + CONFIG::D));
1206  Aas.resize(ptsExtend_.size(), (1 + data_points.size() + CONFIG::D));
1207 
1208  //set Css
1209  // Define intermediate matrix for performance purpose
1210  linalg::sparse_matrix<INT, REAL> Css_coo((1 + data_points.size() + CONFIG::D),(1 + data_points.size() + CONFIG::D),"COO");
1211 
1212  for ( size_t i = 0; i < data_points.size(); i++ ) {
1213  for ( size_t j = i; j < data_points.size(); j++ ) {
1214  auto d = norm(data_points[i].first - data_points[j].first);
1215 
1216  if ( d < r_ ) {
1217  REAL w = rbf(d);
1218  Css_coo.set_value((i + CONFIG::D + 1), (j + CONFIG::D + 1), w, false);
1219 // Css.set_value((i + CONFIG::D + 1), (j + CONFIG::D + 1), w);
1220 
1221  if ( i != j )
1222  Css_coo.set_value((j + CONFIG::D + 1), (i + CONFIG::D + 1), w, false);
1223 // Css.set_value((j + CONFIG::D + 1), (i + CONFIG::D + 1), w);
1224  }
1225  }
1226  }
1227  Css = Css_coo;
1228 
1229  //set Aas
1230  // Define intermediate matrix for performance purpose
1231  linalg::sparse_matrix<INT, REAL> Aas_coo(ptsExtend_.size(),(1 + data_points.size() + CONFIG::D),"COO");
1232 
1233  for ( size_t i = 0; i < ptsExtend_.size(); i++ ) {
1234  for ( size_t j = 0; j < data_points.size(); j++ ) {
1235  auto d = norm(ptsExtend_[i] - data_points[j].first);
1236 
1237  if ( d < r_ ) {
1238 
1239  Aas_coo.set_value(i, (j + CONFIG::D + 1), rbf(d), false);
1240 // Aas.set_value(i, (j + CONFIG::D + 1), rbf(d));
1241  }
1242  }
1243  }
1244  Aas = Aas_coo;
1245 
1246  linalg::sparse_matrix<INT,REAL> Aas_trans = Aas.transpose();
1247 
1248  linalg::sparse_matrix<INT, REAL> H_more;
1249 
1250  if (precond_ == 0) {
1251  linalg::conjugate_gradient<INT, REAL> cg(Css, Aas_trans, cgSolveTol_, cgMaxIter_);
1252  iterErrorReturn = cg.solve();
1253  H_more = cg.getSolution();
1254  }
1255  else if (precond_ == 1) {
1256  linalg::diagonal_preconditioner<INT, REAL> M(Css);
1257  linalg::conjugate_gradient<INT, REAL> cg(Css, Aas_trans, cgSolveTol_, cgMaxIter_, &M);
1258  iterErrorReturn = cg.solve();
1259  H_more = cg.getSolution();
1260  }
1261  else {
1262  std::cerr
1263  << "MUI Error [sampler_rbf.h]: Invalid CG Preconditioner function number ("
1264  << precond_ << ")" << std::endl
1265  << "Please set the CG Preconditioner function number (precond_) as: "
1266  << std::endl << "precond_=0 (No Preconditioner); "
1267  << std::endl << "precond_=1 (Diagonal Preconditioner); " << std::endl;
1268  std::abort();
1269  }
1270 
1271  if (DEBUG) {
1272  std::cout << "MUI [sampler_rbf.h]: #iterations of H_more: " << iterErrorReturn.first
1273  << ". Error of H_more: " << iterErrorReturn.second
1274  << std::endl;
1275  }
1276 
1277  errorReturn = iterErrorReturn.second;
1278 
1279  if ( smoothing ) {
1280  for ( size_t i = 0; i < data_points.size(); i++ ) {
1281  for (size_t j = 0; j < ptsExtend_.size(); j++ ) {
1282  H_toSmooth_.set_value(j, i, H_more.get_value((i + CONFIG::D + 1), j));
1283  }
1284  }
1285  }
1286  else {
1287  for ( size_t i = 0; i < data_points.size(); i++ ) {
1288  for ( size_t j = 0; j < ptsExtend_.size(); j++ ) {
1289  H_.set_value(j, i, H_more.get_value((i + CONFIG::D + 1), j));
1290  }
1291  }
1292  }
1293 
1294  if ( smoothing ) {
1295  for ( size_t row = 0; row < ptsExtend_.size(); row++ ) {
1296  for ( size_t j = 0; j < data_points.size(); j++ ) {
1297  REAL h_j_sum = 0.;
1298  REAL f_sum = 0.;
1299  for ( size_t k = 0; k < MP; k++ ) {
1300  INT row_k = connectivityAA_[row][k];
1301  if ( row_k == static_cast<INT>(row) )
1302  std::cerr << "Invalid row_k value: " << row_k << std::endl;
1303  else
1304  h_j_sum += std::pow(dist_h_i(row, row_k), -2.);
1305  }
1306 
1307  for ( size_t k = 0; k < MP; k++ ) {
1308  INT row_k = connectivityAA_[row][k];
1309  if ( row_k == static_cast<INT>(row) )
1310  std::cerr << "Invalid row_k value: " << row_k << std::endl;
1311  else {
1312  REAL w_i = ((std::pow(dist_h_i(row, row_k), -2.)) / (h_j_sum));
1313  f_sum += w_i * H_toSmooth_.get_value(row_k, j);
1314  }
1315  }
1316  H_.set_value(row, j, 0.5 * (f_sum + H_toSmooth_.get_value(row, j)));
1317  }
1318  }
1319  }
1320  }
1321 
1322  return errorReturn;
1323  }
1324 
1325  template<template<typename, typename > class CONTAINER>
1326  inline REAL buildMatrixConservative(
1327  const CONTAINER<ITYPE, CONFIG> &data_points, const size_t NP,
1328  const size_t MP, bool smoothing, bool pou) const {
1329  REAL errorReturn = 0;
1330  std::pair<INT, REAL> iterErrorReturn(0, 0);
1331  if( pou ) { // Using partitioned approach
1332  for (size_t row = 0; row < data_points.size(); row++) {
1333  linalg::sparse_matrix<INT, REAL> Css; //< Matrix of radial basis function evaluations between prescribed points
1334  linalg::sparse_matrix<INT, REAL> Aas; //< Matrix of RBF evaluations between prescribed and interpolation points
1335 
1336  Css.resize((1 + NP + CONFIG::D), (1 + NP + CONFIG::D));
1337  Aas.resize((1 + NP + CONFIG::D), 1);
1338 
1339  //set Css
1340  // Define intermediate matrix for performance purpose
1341  linalg::sparse_matrix<INT, REAL> Css_coo((1 + NP + CONFIG::D),(1 + NP + CONFIG::D),"COO");
1342 
1343  for (size_t i = 0; i < NP; i++) {
1344  for (size_t j = i; j < NP; j++) {
1345  INT glob_i = connectivityAB_[row][i];
1346  INT glob_j = connectivityAB_[row][j];
1347 
1348  auto d = norm(ptsExtend_[glob_i] - ptsExtend_[glob_j]);
1349 
1350  if (d < r_) {
1351  REAL w = rbf(d);
1352  Css_coo.set_value(i, j, w, false);
1353 // Css.set_value(i, j, w);
1354  if (i != j) {
1355  Css_coo.set_value(j, i, w, false);
1356 // Css.set_value(j, i, w);
1357  }
1358  }
1359  }
1360  }
1361 
1362  for (size_t i = 0; i < NP; i++) {
1363  Css_coo.set_value(i, NP, 1, false);
1364 // Css.set_value(i, NP, 1);
1365  Css_coo.set_value(NP, i, 1, false);
1366 // Css.set_value(NP, i, 1);
1367 
1368  INT glob_i = connectivityAB_[row][i];
1369 
1370  for (INT dim = 0; dim < CONFIG::D; dim++) {
1371  Css_coo.set_value(i, (NP + dim + 1), ptsExtend_[glob_i][dim], false);
1372 // Css.set_value(i, (NP + dim + 1), ptsExtend_[glob_i][dim]);
1373  Css_coo.set_value((NP + dim + 1), i, ptsExtend_[glob_i][dim], false);
1374 // Css.set_value((NP + dim + 1), i, ptsExtend_[glob_i][dim]);
1375  }
1376  }
1377 
1378  Css = Css_coo;
1379 
1380  //set Aas
1381  // Define intermediate matrix for performance purpose
1382  linalg::sparse_matrix<INT, REAL> Aas_coo((1 + NP + CONFIG::D),1,"COO");
1383 
1384  for (size_t j = 0; j < NP; j++) {
1385  INT glob_j = connectivityAB_[row][j];
1386 
1387  auto d = norm(data_points[row].first - ptsExtend_[glob_j]);
1388 
1389  if (d < r_) {
1390  Aas_coo.set_value(j, 0, rbf(d), false);
1391 // Aas.set_value(j, 0, rbf(d));
1392  }
1393  }
1394 
1395  Aas_coo.set_value(NP, 0, 1, false);
1396 // Aas.set_value(NP, 0, 1);
1397 
1398  for (int dim = 0; dim < CONFIG::D; dim++) {
1399  Aas_coo.set_value((NP + dim + 1), 0, data_points[row].first[dim], false);
1400 // Aas.set_value((NP + dim + 1), 0, data_points[row].first[dim]);
1401  }
1402  Aas = Aas_coo;
1403 
1404  linalg::sparse_matrix<INT, REAL> H_i;
1405 
1406  if (precond_ == 0) {
1407  linalg::conjugate_gradient<INT, REAL> cg(Css, Aas, cgSolveTol_, cgMaxIter_);
1408  iterErrorReturn = cg.solve();
1409  H_i = cg.getSolution();
1410  }
1411  else if (precond_ == 1) {
1412  linalg::diagonal_preconditioner<INT, REAL> M(Css);
1413  linalg::conjugate_gradient<INT, REAL> cg(Css, Aas, cgSolveTol_, cgMaxIter_, &M);
1414  iterErrorReturn = cg.solve();
1415  H_i = cg.getSolution();
1416  }
1417  else {
1418  std::cerr
1419  << "MUI Error [sampler_rbf.h]: Invalid CG Preconditioner function number ("
1420  << precond_ << ")" << std::endl
1421  << "Please set the CG Preconditioner function number (precond_) as: "
1422  << std::endl << "precond_=0 (No Preconditioner); "
1423  << std::endl << "precond_=1 (Diagonal Preconditioner); " << std::endl;
1424  std::abort();
1425  }
1426 
1427  if (DEBUG) {
1428  std::cout << "MUI [sampler_rbf.h]: #iterations of H_i: " << iterErrorReturn.first
1429  << ". Error of H_i: " << iterErrorReturn.second
1430  << std::endl;
1431  }
1432 
1433  errorReturn += iterErrorReturn.second;
1434 
1435  if (smoothing) {
1436  for (size_t j = 0; j < NP; j++) {
1437  INT glob_j = connectivityAB_[row][j];
1438  H_toSmooth_.set_value(glob_j, row, H_i.get_value(j, 0));
1439  }
1440  }
1441  else {
1442  for (size_t j = 0; j < NP; j++) {
1443  INT glob_j = connectivityAB_[row][j];
1444  H_.set_value(glob_j, row, H_i.get_value(j, 0));
1445  }
1446  }
1447  }
1448 
1449  errorReturn /= static_cast<REAL>(data_points.size());
1450 
1451  if (smoothing) {
1452  for (size_t row = 0; row < data_points.size(); row++) {
1453  for (size_t j = 0; j < NP; j++) {
1454  INT row_i = connectivityAB_[row][j];
1455  REAL h_j_sum = 0.;
1456  REAL f_sum = 0.;
1457 
1458  for (size_t k = 0; k < MP; k++) {
1459  INT row_k = connectivityAA_[row_i][k];
1460  if (row_k == static_cast<INT>(row_i)) {
1461  std::cerr << "MUI Error [sampler_rbf.h]: Invalid row_k value: " << row_k
1462  << std::endl;
1463  }
1464  else
1465  h_j_sum += std::pow(dist_h_i(row_i, row_k), -2.);
1466  }
1467 
1468  for (size_t k = 0; k < MP; k++) {
1469  INT row_k = connectivityAA_[row_i][k];
1470  if (row_k == static_cast<INT>(row_i)) {
1471  std::cerr << "MUI Error [sampler_rbf.h]: Invalid row_k value: " << row_k
1472  << std::endl;
1473  }
1474  else {
1475  REAL w_i = ((std::pow(dist_h_i(row_i, row_k), -2.))
1476  / (h_j_sum));
1477  f_sum += w_i * H_toSmooth_.get_value(row_k, row);
1478  }
1479  }
1480  H_.set_value(row_i, row,
1481  (0.5 * (f_sum + H_toSmooth_.get_value(row_i, row))));
1482  }
1483  }
1484  }
1485  }
1486  else { // Not using partitioned approach
1487  linalg::sparse_matrix<INT,REAL> Css; //< Matrix of radial basis function evaluations between prescribed points
1488  linalg::sparse_matrix<INT,REAL> Aas; //< Matrix of RBF evaluations between prescribed and interpolation points
1489 
1490  Css.resize((1 + ptsExtend_.size() + CONFIG::D), (1 + ptsExtend_.size() + CONFIG::D));
1491  Aas.resize(data_points.size(), (1 + ptsExtend_.size() + CONFIG::D));
1492 
1493  //set Css
1494  // Define intermediate vectors for performance purpose
1495  linalg::sparse_matrix<INT, REAL> Css_coo((1 + ptsExtend_.size() + CONFIG::D),(1 + ptsExtend_.size() + CONFIG::D),"COO");
1496 
1497  for ( size_t i = 0; i < ptsExtend_.size(); i++ ) {
1498  for ( size_t j = i; j < ptsExtend_.size(); j++ ) {
1499  auto d = norm(ptsExtend_[i] - ptsExtend_[j]);
1500 
1501  if ( d < r_ ) {
1502  REAL w = rbf(d);
1503  Css_coo.set_value((i + CONFIG::D + 1), (j + CONFIG::D + 1), w, false);
1504 // Css.set_value((i + CONFIG::D + 1), (j + CONFIG::D + 1), w);
1505 
1506  if ( i != j ) {
1507  Css_coo.set_value((j + CONFIG::D + 1), (i + CONFIG::D + 1), w, false);
1508 // Css.set_value((j + CONFIG::D + 1), (i + CONFIG::D + 1), w);
1509  }
1510  }
1511  }
1512  }
1513  Css = Css_coo;
1514 
1515  //set Aas
1516  // Define intermediate vectors for performance purpose
1517  linalg::sparse_matrix<INT, REAL> Aas_coo(data_points.size(),(1 + ptsExtend_.size() + CONFIG::D),"COO");
1518 
1519  for ( size_t i = 0; i < data_points.size(); i++ ) {
1520  for ( size_t j = 0; j < ptsExtend_.size(); j++ ) {
1521  auto d = norm(data_points[i].first - ptsExtend_[j]);
1522 
1523  if ( d < r_ ) {
1524  Aas_coo.set_value(i, (j + CONFIG::D + 1), rbf(d), false);
1525 // Aas.set_value(i, (j + CONFIG::D + 1), rbf(d));
1526  }
1527  }
1528  }
1529  Aas = Aas_coo;
1530 
1531  linalg::sparse_matrix<INT,REAL> Aas_trans = Aas.transpose();
1532 
1533  linalg::sparse_matrix<INT, REAL> H_more;
1534 
1535  if (precond_ == 0) {
1536  linalg::conjugate_gradient<INT, REAL> cg(Css, Aas_trans, cgSolveTol_, cgMaxIter_);
1537  iterErrorReturn = cg.solve();
1538  H_more = cg.getSolution();
1539  }
1540  else if (precond_ == 1) {
1541  linalg::diagonal_preconditioner<INT, REAL> M(Css);
1542  linalg::conjugate_gradient<INT, REAL> cg(Css, Aas_trans, cgSolveTol_, cgMaxIter_, &M);
1543  iterErrorReturn = cg.solve();
1544  H_more = cg.getSolution();
1545  }
1546  else {
1547  std::cerr
1548  << "MUI Error [sampler_rbf.h]: Invalid CG Preconditioner function number ("
1549  << precond_ << ")" << std::endl
1550  << "Please set the CG Preconditioner function number (precond_) as: "
1551  << std::endl << "precond_=0 (No Preconditioner); "
1552  << std::endl << "precond_=1 (Diagonal Preconditioner); " << std::endl;
1553  std::abort();
1554  }
1555 
1556  if (DEBUG) {
1557  std::cout << "MUI [sampler_rbf.h]: #iterations of H_more: " << iterErrorReturn.first
1558  << ". Error of H_more: " << iterErrorReturn.second
1559  << std::endl;
1560  }
1561 
1562  errorReturn = iterErrorReturn.second;
1563 
1564  if ( smoothing ) {
1565  for ( size_t i = 0; i < ptsExtend_.size(); i++ ) {
1566  for (size_t j = 0; j < data_points.size(); j++ ) {
1567  H_toSmooth_.set_value(i, j, H_more.get_value((i + CONFIG::D + 1), j));
1568  }
1569  }
1570  }
1571  else {
1572  for ( size_t i = 0; i < ptsExtend_.size(); i++ ) {
1573  for ( size_t j = 0; j < data_points.size(); j++ ) {
1574  H_.set_value(i, j, H_more.get_value((i + CONFIG::D + 1), j));
1575  }
1576  }
1577  }
1578 
1579  if ( smoothing ) {
1580  for ( size_t row = 0; row < ptsExtend_.size(); row++ ) {
1581  for ( size_t j = 0; j < data_points.size(); j++ ) {
1582  REAL h_j_sum = 0.;
1583  REAL f_sum = 0.;
1584  for ( size_t k = 0; k < MP; k++ ) {
1585  INT row_k = connectivityAA_[row][k];
1586  if ( row_k == static_cast<INT>(row) )
1587  std::cerr << "Invalid row_k value: " << row_k << std::endl;
1588  else
1589  h_j_sum += std::pow(dist_h_i(row, row_k), -2.);
1590  }
1591 
1592  for ( size_t k = 0; k < MP; k++ ) {
1593  INT row_k = connectivityAA_[row][k];
1594  if ( row_k == static_cast<INT>(row) )
1595  std::cerr << "Invalid row_k value: " << row_k << std::endl;
1596  else {
1597  REAL w_i = ((std::pow(dist_h_i(row, row_k), -2.)) / (h_j_sum));
1598  f_sum += w_i * H_toSmooth_.get_value(row_k, j);
1599  }
1600  }
1601  H_.set_value(row, j, 0.5 * (f_sum + H_toSmooth_.get_value(row, j)));
1602  }
1603  }
1604  }
1605  }
1606 
1607  return errorReturn;
1608  }
1609 
1610  template<template<typename, typename > class CONTAINER>
1611  inline void buildConnectivityConsistent(const CONTAINER<ITYPE, CONFIG> &data_points, const size_t NP, bool writeMatrix,
1612  const std::string& fileAddress) const {
1613  std::ofstream outputFileCAB;
1614  std::ofstream outputFileRemotePoints;
1615  if (writeMatrix) {
1616  outputFileCAB.open(fileAddress + "/connectivityAB.dat");
1617 
1618  if (!outputFileCAB) {
1619  std::cerr
1620  << "MUI Error [sampler_rbf.h]: Could not locate the file address on the connectivityAB.dat"
1621  << std::endl;
1622  }
1623  else {
1624  outputFileCAB
1625  << "// ************************************************************************************************";
1626  outputFileCAB << "\n";
1627  outputFileCAB
1628  << "// **** This is the 'connectivityAB.dat' file of the RBF spatial sampler of the MUI library";
1629  outputFileCAB << "\n";
1630  outputFileCAB
1631  << "// **** This file contains the entire matrix of the Point Connectivity Matrix (N).";
1632  outputFileCAB << "\n";
1633  outputFileCAB
1634  << "// **** The file uses the Comma-Separated Values (CSV) format with ASCII for the entire N matrix";
1635  outputFileCAB << "\n";
1636  outputFileCAB
1637  << "// ************************************************************************************************";
1638  outputFileCAB << "\n";
1639  outputFileCAB << "// ";
1640  outputFileCAB << "\n";
1641  }
1642 
1643  outputFileRemotePoints.open(fileAddress + "/remotePoints.dat");
1644 
1645  if (!outputFileRemotePoints) {
1646  std::cerr
1647  << "MUI Error [sampler_rbf.h]: Could not locate the file address on the remotePoints.dat"
1648  << std::endl;
1649  }
1650  else {
1651  outputFileRemotePoints
1652  << "// ************************************************************************************************";
1653  outputFileRemotePoints << "\n";
1654  outputFileRemotePoints
1655  << "// **** This is the 'remotePoints.dat' file of the RBF spatial sampler of the MUI library";
1656  outputFileRemotePoints << "\n";
1657  outputFileRemotePoints
1658  << "// **** This file contains the remote points in the correct order that build the coupling matrix H.";
1659  outputFileRemotePoints << "\n";
1660  outputFileRemotePoints
1661  << "// **** The file uses the Comma-Separated Values (CSV) format with ASCII for the entire N matrix";
1662  outputFileRemotePoints << "\n";
1663  outputFileRemotePoints
1664  << "// ************************************************************************************************";
1665  outputFileRemotePoints << "\n";
1666  outputFileRemotePoints << "// ";
1667  outputFileRemotePoints << "\n";
1668  }
1669  }
1670 
1671  connectivityAB_.resize(ptsExtend_.size());
1672 
1673  for (size_t i = 0; i < ptsExtend_.size(); i++) {
1674  INT pointsCount = 0;
1675  for (size_t n = 0; n < NP; n++) {
1676  REAL cur = std::numeric_limits<REAL>::max();
1677  INT bestj = -1;
1678  for (size_t j = 0; j < data_points.size(); j++) {
1679  auto added = std::find_if(connectivityAB_[i].begin(),
1680  connectivityAB_[i].end(), [j](INT k) {
1681  return static_cast<size_t>(k) == j;
1682  });
1683 
1684  if (added != connectivityAB_[i].end())
1685  continue;
1686 
1687  auto d = normsq(ptsExtend_[i] - data_points[j].first);
1688  if (d < cur) {
1689  cur = d;
1690  bestj = j;
1691  }
1692 
1693  if (n == 0 && d < twor_)
1694  pointsCount++;
1695  }
1696 
1697  connectivityAB_[i].emplace_back(bestj);
1698 
1699  if (writeMatrix && (n < NP - 1))
1700  outputFileCAB << bestj << ",";
1701  else if (writeMatrix)
1702  outputFileCAB << bestj;
1703  }
1704 
1705  if (writeMatrix && i < ptsExtend_.size() - 1)
1706  outputFileCAB << '\n';
1707  }
1708 
1709  for (size_t i = 0; i < data_points.size(); i++) {
1710 
1711  point_type pointTemp;
1712  for (INT dim = 0; dim < CONFIG::D; dim++) {
1713  pointTemp[dim] = data_points[i].first[dim];
1714  }
1715 
1716  remote_pts_.push_back(pointTemp);
1717 
1718  if (writeMatrix) {
1719  for (INT dim = 0; dim < CONFIG::D; dim++) {
1720  outputFileRemotePoints << remote_pts_[i][dim];
1721  if (dim == (CONFIG::D - 1)) {
1722  outputFileRemotePoints << '\n';
1723  } else {
1724  outputFileRemotePoints << ",";
1725  }
1726  }
1727  }
1728  }
1729 
1730  if (writeMatrix) {
1731  outputFileCAB.close();
1732  outputFileRemotePoints.close();
1733  }
1734  }
1735 
1736  template<template<typename, typename > class CONTAINER>
1737  inline void buildConnectivityConservative(const CONTAINER<ITYPE, CONFIG> &data_points, const size_t NP, bool writeMatrix,
1738  const std::string& fileAddress) const {
1739  std::ofstream outputFileCAB;
1740  std::ofstream outputFileRemotePoints;
1741  if (writeMatrix) {
1742  outputFileCAB.open(fileAddress + "/connectivityAB.dat");
1743 
1744  if (!outputFileCAB) {
1745  std::cerr
1746  << "MUI Error [sampler_rbf.h]: Could not locate the file address on the connectivityAB.dat"
1747  << std::endl;
1748  }
1749  else {
1750  outputFileCAB
1751  << "// ************************************************************************************************";
1752  outputFileCAB << "\n";
1753  outputFileCAB
1754  << "// **** This is the 'connectivityAB.dat' file of the RBF spatial sampler of the MUI library";
1755  outputFileCAB << "\n";
1756  outputFileCAB
1757  << "// **** This file contains the entire matrix of the Point Connectivity Matrix (N).";
1758  outputFileCAB << "\n";
1759  outputFileCAB
1760  << "// **** The file uses the Comma-Separated Values (CSV) format with ASCII for the entire N matrix";
1761  outputFileCAB << "\n";
1762  outputFileCAB
1763  << "// ************************************************************************************************";
1764  outputFileCAB << "\n";
1765  outputFileCAB << "// ";
1766  outputFileCAB << "\n";
1767  }
1768 
1769  outputFileRemotePoints.open(fileAddress + "/remotePoints.dat");
1770 
1771  if (!outputFileRemotePoints) {
1772  std::cerr
1773  << "MUI Error [sampler_rbf.h]: Could not locate the file address on the remotePoints.dat"
1774  << std::endl;
1775  }
1776  else {
1777  outputFileRemotePoints
1778  << "// ************************************************************************************************";
1779  outputFileRemotePoints << "\n";
1780  outputFileRemotePoints
1781  << "// **** This is the 'remotePoints.dat' file of the RBF spatial sampler of the MUI library";
1782  outputFileRemotePoints << "\n";
1783  outputFileRemotePoints
1784  << "// **** This file contains the remote points in the correct order that build the coupling matrix H.";
1785  outputFileRemotePoints << "\n";
1786  outputFileRemotePoints
1787  << "// **** The file uses the Comma-Separated Values (CSV) format with ASCII for the entire N matrix";
1788  outputFileRemotePoints << "\n";
1789  outputFileRemotePoints
1790  << "// ************************************************************************************************";
1791  outputFileRemotePoints << "\n";
1792  outputFileRemotePoints << "// ";
1793  outputFileRemotePoints << "\n";
1794  }
1795  }
1796 
1797  connectivityAB_.resize(data_points.size());
1798 
1799  for (size_t i = 0; i < data_points.size(); i++) {
1800  INT pointsCount = 0;
1801  for (size_t n = 0; n < NP; n++) {
1802  REAL cur = std::numeric_limits<REAL>::max();
1803  INT bestj = -1;
1804  for (size_t j = 0; j < ptsExtend_.size(); j++) {
1805  auto added = std::find_if(connectivityAB_[i].begin(),
1806  connectivityAB_[i].end(), [j](INT k) {
1807  return static_cast<size_t>(k) == j;
1808  });
1809 
1810  if (added != connectivityAB_[i].end())
1811  continue;
1812 
1813  auto d = normsq(data_points[i].first - ptsExtend_[j]);
1814  if (d < cur) {
1815  cur = d;
1816  bestj = j;
1817  }
1818 
1819  if (n == 0 && d < twor_)
1820  pointsCount++;
1821  }
1822 
1823  connectivityAB_[i].emplace_back(bestj);
1824 
1825  if (writeMatrix && n < NP - 1)
1826  outputFileCAB << bestj << ",";
1827  else if (writeMatrix)
1828  outputFileCAB << bestj;
1829  }
1830 
1831  if (writeMatrix && i < ptsExtend_.size() - 1)
1832  outputFileCAB << '\n';
1833 
1834  point_type pointTemp;
1835  for (INT dim = 0; dim < CONFIG::D; dim++) {
1836  pointTemp[dim] = data_points[i].first[dim];
1837  }
1838  remote_pts_.push_back(pointTemp);
1839 
1840  if (writeMatrix) {
1841  for (INT dim = 0; dim < CONFIG::D; dim++) {
1842  outputFileRemotePoints << remote_pts_[i][dim];
1843  if (dim == (CONFIG::D - 1)) {
1844  outputFileRemotePoints << '\n';
1845  } else {
1846  outputFileRemotePoints << ",";
1847  }
1848  }
1849  }
1850  }
1851 
1852  if (writeMatrix) {
1853  outputFileCAB.close();
1854  outputFileRemotePoints.close();
1855  }
1856  }
1857 
1858  void buildConnectivitySmooth(const size_t MP, bool writeMatrix, const std::string& fileAddress) const {
1859  std::ofstream outputFileCAA;
1860 
1861  if (writeMatrix) {
1862  outputFileCAA.open(fileAddress + "/connectivityAA.dat");
1863 
1864  if (!outputFileCAA) {
1865  std::cerr
1866  << "MUI Error [sampler_rbf.h]: Could not locate the file address on the connectivityAA.dat!"
1867  << std::endl;
1868  }
1869  else {
1870  outputFileCAA
1871  << "// ************************************************************************************************";
1872  outputFileCAA << "\n";
1873  outputFileCAA
1874  << "// **** This is the 'connectivityAA.dat' file of the RBF spatial sampler of the MUI library";
1875  outputFileCAA << "\n";
1876  outputFileCAA
1877  << "// **** This file contains the entire matrix of the Point Connectivity Matrix (M) (for smoothing).";
1878  outputFileCAA << "\n";
1879  outputFileCAA
1880  << "// **** The file uses the Comma-Separated Values (CSV) format with ASCII for the entire N matrix";
1881  outputFileCAA << "\n";
1882  outputFileCAA
1883  << "// ************************************************************************************************";
1884  outputFileCAA << "\n";
1885  outputFileCAA << "// ";
1886  outputFileCAA << "\n";
1887  }
1888  }
1889 
1890  connectivityAA_.resize(ptsExtend_.size());
1891 
1892  for (size_t i = 0; i < ptsExtend_.size(); i++) {
1893  for (size_t n = 0; n < MP; n++) {
1894  REAL cur = std::numeric_limits<REAL>::max();
1895  INT bestj = -1;
1896  for (size_t j = 0; j < ptsExtend_.size(); j++) {
1897  if (i == j)
1898  continue;
1899 
1900  auto added = std::find_if(connectivityAA_[i].begin(),
1901  connectivityAA_[i].end(), [j](INT i) {
1902  return static_cast<size_t>(i) == j;
1903  });
1904 
1905  if (added != connectivityAA_[i].end())
1906  continue;
1907 
1908  auto d = normsq(ptsExtend_[i] - ptsExtend_[j]);
1909  if (d < cur) {
1910  cur = d;
1911  bestj = j;
1912  }
1913  }
1914 
1915  connectivityAA_[i].emplace_back(bestj);
1916 
1917  if (writeMatrix && n < MP - 1)
1918  outputFileCAA << bestj << ",";
1919  else if (writeMatrix)
1920  outputFileCAA << bestj;
1921  }
1922 
1923  if (writeMatrix && i < ptsExtend_.size() - 1)
1924  outputFileCAA << '\n';
1925  }
1926 
1927  if (writeMatrix)
1928  outputFileCAA.close();
1929  }
1930 
1931  //Radial basis function for two points
1932  inline REAL rbf(point_type x1, point_type x2) const {
1933  auto d = norm(x1 - x2);
1934  return rbf(d);
1935  }
1936 
1937  //Radial basis function for calculated distance
1938  inline REAL rbf(REAL d) const {
1939  switch (basisFunc_) {
1940  case 0:
1941  //Gaussian
1942  return (d < r_) ? std::exp(-(s_ * s_ * d * d)) : 0.;
1943  case 1:
1944  //Wendland's C0
1945  return std::pow((1. - d), 2.);
1946  case 2:
1947  //Wendland's C2
1948  return (std::pow((1. - d), 4.)) * ((4. * d) + 1.);
1949  case 3:
1950  //Wendland's C4
1951  return (std::pow((1. - d), 6)) * ((35. * d * d) + (18. * d) + 3.);
1952  case 4:
1953  //Wendland's C6
1954  return (std::pow((1. - d), 8.))
1955  * ((32. * d * d * d) + (25. * d * d) + (8. * d) + 1.);
1956  default:
1957  std::cerr
1958  << "MUI Error [sampler_rbf.h]: invalid RBF basis function number ("
1959  << basisFunc_ << ")" << std::endl
1960  << "Please set the RBF basis function number (basisFunc_) as: "
1961  << std::endl << "basisFunc_=0 (Gaussian); " << std::endl
1962  << "basisFunc_=1 (Wendland's C0); " << std::endl
1963  << "basisFunc_=2 (Wendland's C2); " << std::endl
1964  << "basisFunc_=3 (Wendland's C4); " << std::endl
1965  << "basisFunc_=4 (Wendland's C6); " << std::endl;
1966  return 0;
1967  }
1968  }
1969 
1971  inline REAL dist_h_i(INT ptsExtend_i, INT ptsExtend_j) const {
1972  switch (CONFIG::D) {
1973  case 1:
1974  return std::sqrt((std::pow((ptsExtend_[ptsExtend_i][0]
1975  - ptsExtend_[ptsExtend_j][0]), 2.)));
1976  case 2:
1977  return std::sqrt((std::pow((ptsExtend_[ptsExtend_i][0]
1978  - ptsExtend_[ptsExtend_j][0]), 2.))
1979  + (std::pow((ptsExtend_[ptsExtend_i][1]
1980  - ptsExtend_[ptsExtend_j][1]), 2.)));
1981  case 3:
1982  return std::sqrt((std::pow((ptsExtend_[ptsExtend_i][0]
1983  - ptsExtend_[ptsExtend_j][0]), 2.))
1984  + (std::pow((ptsExtend_[ptsExtend_i][1]
1985  - ptsExtend_[ptsExtend_j][1]), 2.))
1986  + (std::pow((ptsExtend_[ptsExtend_i][2]
1987  - ptsExtend_[ptsExtend_j][2]), 2.)));
1988  default:
1989  std::cerr << "MUI Error [sampler_rbf.h]: CONFIG::D must equal 1-3" << std::endl;
1990  return 0.;
1991  }
1992  }
1993 
1994 protected:
1995  REAL r_;
1996  const std::vector<point_type> pts_; //< Local points
1997  const INT basisFunc_;
1998  const bool conservative_;
1999  const bool consistent_;
2000  const bool smoothFunc_;
2001  const bool generateMatrix_;
2002  const std::string writeFileAddress_;
2003  INT precond_;
2004  mutable bool initialised_;
2005  MPI_Comm local_mpi_comm_world_;
2006  mutable INT CABrow_;
2007  mutable INT CABcol_;
2008  mutable INT Hrow_;
2009  mutable INT Hcol_;
2010  bool pouEnabled_;
2011  REAL cgSolveTol_;
2012  INT cgMaxIter_;
2013 
2014  mutable size_t N_sp_;
2015  mutable size_t M_ap_;
2016  int local_rank_;
2017  int local_size_;
2018 
2019  REAL twor_;
2020  REAL s_;
2021  mutable INT CAArow_;
2022  mutable INT CAAcol_;
2023  mutable INT remote_pts_num_;
2024  mutable INT remote_pts_dim_;
2025  mutable std::vector<point_type> ptsGhost_; //< Local ghost points
2026  mutable std::vector<point_type> ptsExtend_; //< Extended local points, i.e. local points and ghost local points
2027  mutable std::vector<std::vector<INT> > connectivityAB_;
2028  mutable std::vector<std::vector<INT> > connectivityAA_;
2029  mutable linalg::sparse_matrix<INT, REAL> H_; //< Transformation Matrix
2030  mutable linalg::sparse_matrix<INT, REAL> H_toSmooth_;
2031  mutable std::vector<point_type> remote_pts_; //< Remote points
2032 };
2033 } // mui
2034 
2035 #endif /* MUI_SAMPLER_RBF_H_ */
mui::geometry::any_shape
Definition: geometry.h:92
mui::linalg::sparse_matrix< INT, REAL >
mui::linalg::sparse_matrix::get_rows
ITYPE get_rows() const
Definition: matrix_io_info.h:579
mui::linalg::sparse_matrix::set_zero
void set_zero()
Definition: matrix_manipulation.h:418
mui::linalg::sparse_matrix::set_value
void set_value(ITYPE, ITYPE, VTYPE, bool=true)
Definition: matrix_manipulation.h:292
mui::linalg::sparse_matrix::get_value
VTYPE get_value(ITYPE, ITYPE) const
Definition: matrix_io_info.h:523
mui::linalg::sparse_matrix::resize
void resize(ITYPE, ITYPE)
Definition: matrix_manipulation.h:62
mui::linalg::sparse_matrix::get_cols
ITYPE get_cols() const
Definition: matrix_io_info.h:585
mui::sampler_rbf
Definition: sampler_rbf.h:63
mui::sampler_rbf::initialised_
bool initialised_
Definition: sampler_rbf.h:2004
mui::sampler_rbf::cgSolveTol_
REAL cgSolveTol_
Definition: sampler_rbf.h:2011
mui::sampler_rbf::ITYPE
I_TP ITYPE
Definition: sampler_rbf.h:67
mui::sampler_rbf::r_
REAL r_
Definition: sampler_rbf.h:1995
mui::sampler_rbf::connectivityAA_
std::vector< std::vector< INT > > connectivityAA_
Definition: sampler_rbf.h:2028
mui::sampler_rbf::getGhostPointsToSend
std::pair< std::vector< std::pair< INT, INT > >, std::vector< std::pair< INT, std::vector< point_type > > > > getGhostPointsToSend(std::pair< point_type, point_type > lbbExtend, MPI_Comm local_world, int local_rank, int local_size) const
Definition: sampler_rbf.h:565
mui::sampler_rbf::N_sp_
size_t N_sp_
Definition: sampler_rbf.h:2014
mui::sampler_rbf::addFetchPoint
void addFetchPoint(point_type pt)
Definition: sampler_rbf.h:308
mui::sampler_rbf::CABcol_
INT CABcol_
Definition: sampler_rbf.h:2007
mui::sampler_rbf::basisFunc_
const INT basisFunc_
Definition: sampler_rbf.h:1997
mui::sampler_rbf::preSetFetchPointsExtend
void preSetFetchPointsExtend(std::vector< point_type > &pts)
Definition: sampler_rbf.h:303
mui::sampler_rbf::EXCEPTION
typename CONFIG::EXCEPTION EXCEPTION
Definition: sampler_rbf.h:71
mui::sampler_rbf::connectivityAB_
std::vector< std::vector< INT > > connectivityAB_
Definition: sampler_rbf.h:2027
mui::sampler_rbf::addFetchPointGhost
void addFetchPointGhost(point_type pt)
Definition: sampler_rbf.h:318
mui::sampler_rbf::localExtendBoundingBox
std::pair< point_type, point_type > localExtendBoundingBox(std::pair< point_type, point_type > lbb, REAL r) const
Definition: sampler_rbf.h:553
mui::sampler_rbf::distributeGhostPoints
std::vector< point_type > distributeGhostPoints(std::vector< std::pair< INT, INT >> ghostPointsCountToSend, std::vector< std::pair< INT, std::vector< point_type >>> ghostPointsToSend, MPI_Comm local_world, int local_rank) const
Definition: sampler_rbf.h:748
mui::sampler_rbf::REAL
typename CONFIG::REAL REAL
Definition: sampler_rbf.h:68
mui::sampler_rbf::CAAcol_
INT CAAcol_
Definition: sampler_rbf.h:2022
mui::sampler_rbf::Hrow_
INT Hrow_
Definition: sampler_rbf.h:2008
mui::sampler_rbf::local_size_
int local_size_
Definition: sampler_rbf.h:2017
mui::sampler_rbf::writeFileAddress_
const std::string writeFileAddress_
Definition: sampler_rbf.h:2002
mui::sampler_rbf::support
geometry::any_shape< CONFIG > support(point_type focus, REAL domain_mag) const
Definition: sampler_rbf.h:294
mui::sampler_rbf::pouEnabled_
bool pouEnabled_
Definition: sampler_rbf.h:2010
mui::sampler_rbf::ptsGhost_
std::vector< point_type > ptsGhost_
Definition: sampler_rbf.h:2025
mui::sampler_rbf::preSetFetchPoints
void preSetFetchPoints(std::vector< point_type > &pts)
Definition: sampler_rbf.h:298
mui::sampler_rbf::DEBUG
static const bool DEBUG
Definition: sampler_rbf.h:74
mui::sampler_rbf::CABrow_
INT CABrow_
Definition: sampler_rbf.h:2006
mui::sampler_rbf::consistent_
const bool consistent_
Definition: sampler_rbf.h:1999
mui::sampler_rbf::twor_
REAL twor_
Definition: sampler_rbf.h:2019
mui::sampler_rbf::conservative_
const bool conservative_
Definition: sampler_rbf.h:1998
mui::sampler_rbf::ptsExtend_
std::vector< point_type > ptsExtend_
Definition: sampler_rbf.h:2026
mui::sampler_rbf::remote_pts_num_
INT remote_pts_num_
Definition: sampler_rbf.h:2023
mui::sampler_rbf::cgMaxIter_
INT cgMaxIter_
Definition: sampler_rbf.h:2012
mui::sampler_rbf::point_type
typename CONFIG::point_type point_type
Definition: sampler_rbf.h:70
mui::sampler_rbf::readRBFMatrix
void readRBFMatrix(const std::string &readFileAddress) const
Definition: sampler_rbf.h:323
mui::sampler_rbf::INT
typename CONFIG::INT INT
Definition: sampler_rbf.h:69
mui::sampler_rbf::M_ap_
size_t M_ap_
Definition: sampler_rbf.h:2015
mui::sampler_rbf::addFetchPointExtend
void addFetchPointExtend(point_type pt)
Definition: sampler_rbf.h:313
mui::sampler_rbf::H_
linalg::sparse_matrix< INT, REAL > H_
Definition: sampler_rbf.h:2029
mui::sampler_rbf::localBoundingBox
std::pair< point_type, point_type > localBoundingBox(const std::vector< point_type > pt) const
Definition: sampler_rbf.h:478
mui::sampler_rbf::Hcol_
INT Hcol_
Definition: sampler_rbf.h:2009
mui::sampler_rbf::s_
REAL s_
Definition: sampler_rbf.h:2020
mui::sampler_rbf::precond_
INT precond_
Definition: sampler_rbf.h:2003
mui::sampler_rbf::OTYPE
O_TP OTYPE
Definition: sampler_rbf.h:66
mui::sampler_rbf::filter
OTYPE filter(point_type focus, const CONTAINER< ITYPE, CONFIG > &data_points) const
Definition: sampler_rbf.h:184
mui::sampler_rbf::CAArow_
INT CAArow_
Definition: sampler_rbf.h:2021
mui::sampler_rbf::H_toSmooth_
linalg::sparse_matrix< INT, REAL > H_toSmooth_
Definition: sampler_rbf.h:2030
mui::sampler_rbf::local_mpi_comm_world_
MPI_Comm local_mpi_comm_world_
Definition: sampler_rbf.h:2005
mui::sampler_rbf::local_rank_
int local_rank_
Definition: sampler_rbf.h:2016
mui::sampler_rbf::generateMatrix_
const bool generateMatrix_
Definition: sampler_rbf.h:2001
mui::sampler_rbf::remote_pts_
std::vector< point_type > remote_pts_
Definition: sampler_rbf.h:2031
mui::sampler_rbf::remote_pts_dim_
INT remote_pts_dim_
Definition: sampler_rbf.h:2024
mui::sampler_rbf::smoothFunc_
const bool smoothFunc_
Definition: sampler_rbf.h:2000
mui::sampler_rbf::pts_
const std::vector< point_type > pts_
Definition: sampler_rbf.h:1996
mui::sampler_rbf::sampler_rbf
sampler_rbf(REAL r, const std::vector< point_type > &pts, INT basisFunc=0, bool conservative=false, bool smoothFunc=false, bool generateMatrix=true, const std::string &writeFileAddress=std::string(), REAL cutOff=1e-9, REAL cgSolveTol=1e-6, INT cgMaxIter=0, INT pouSize=50, INT precond=1, MPI_Comm local_comm=MPI_COMM_NULL)
Definition: sampler_rbf.h:126
mui::sampler_rbf::QUIET
static const bool QUIET
Definition: sampler_rbf.h:73
mui::sampler_rbf::facilitateGhostPoints
void facilitateGhostPoints() const
Definition: sampler_rbf.h:848
mui::dim::electrical::e
e
Definition: dim.h:347
mui
Definition: comm.h:54
mui::sum
SCALAR sum(vexpr< E, SCALAR, D > const &u)
Definition: point.h:362
mui::max
SCALAR max(vexpr< E, SCALAR, D > const &u)
Definition: point.h:350
mui::normsq
SCALAR normsq(vexpr< E, SCALAR, D > const &u)
Definition: point.h:380
mui::norm
SCALAR norm(vexpr< E, SCALAR, D > const &u)
Definition: point.h:386