NeoPZ
tpzdohrsubstruct.cpp
Go to the documentation of this file.
1 
6 #include "tpzdohrsubstruct.h"
7 #include <iostream>
8 #include "pzlog.h"
9 
10 #ifdef LOG4CXX
11 static LoggerPtr logger(Logger::getLogger("substruct.dohrsubstruct"));
12 #endif
13 
14 using namespace std;
15 
16 template<class TVar>
17 typename TPZDohrSubstruct<TVar>::EWeightType TPZDohrSubstruct<TVar>::fWeightType = TPZDohrSubstruct<TVar>::CorrectWeight;
18 
19 
20 template<class TVar>
21 TPZDohrSubstruct<TVar>::TPZDohrSubstruct()
22 {
23  //Inicializacao
24 }
25 
26 template<class TVar>
27 TPZDohrSubstruct<TVar>::~TPZDohrSubstruct()
28 {
29  //Limpesa
30 }
31 
32 template<class TVar>
33 void TPZDohrSubstruct<TVar>::ContributeTestV1(TPZFMatrix<TVar> &testV1, int NumCoarse) {
34  /* temp1 will be equal to W(i)*Phi(i) */
35  TPZFMatrix<TVar> temp1(fNEquations,fCoarseIndex.NElements());
36  int i,j;
37  for(i=0;i<fPhiC.Rows();i++) {
38  for(j=0;j<fPhiC.Cols();j++) {
39  temp1(i,j) = fPhiC(i,j)*fWeights[i];
40  }
41  }
42  /* And now temp2 will be equal to temp1*R(ci) */
43  TPZFMatrix<TVar> col(fNEquations,1);
44  TPZFMatrix<TVar> temp2(fNEquations,NumCoarse,0.);
45  for(i=0;i<fCoarseIndex.NElements();i++) {
46  temp1.GetSub(0,i,temp1.Rows(),1,col);
47  temp2.PutSub(0,fCoarseIndex[i],col);
48  }
49  /* And testV1 will be equal to R(i)_trans*temp2 */
50  col.Resize(1,NumCoarse);
51  for(i=0;i<fNEquations;i++) {
52  temp2.GetSub(i,0,1,temp2.Cols(),col);
53  //testV1.PutSub(fGlobalIndex[i],0,col);
54  for(j=0;j<NumCoarse;j++) {
55  testV1(fGlobalIndex[i],j) += col(0,j);
56  }
57  }
58 }
59 
60 template<class TVar>
61 void TPZDohrSubstruct<TVar>::ContributeResidual(TPZFMatrix<TVar> &u, TPZFMatrix<TVar> &r){
62  TPZFMatrix<TVar> ulocal(fNEquations,1,0.);
63  TPZFMatrix<TVar> reslocal(fNEquations,1);
64  int i;
65  int neqs = fGlobalEqs.NElements();
66  for (i=0;i<neqs;i++)
67  {
68  std::pair<int,int> ind = fGlobalEqs[i];
69  ulocal(ind.first,0) = u(ind.second,0);
70  }
71  fStiffness->Residual(ulocal, fLocalLoad, reslocal);
72  for (i=0;i<neqs;i++) {
73  std::pair<int,int> ind = fGlobalEqs[i];
74  r(ind.second,0) += reslocal(ind.first,0);
75  }
76 }
77 
78 template<class TVar>
79 void TPZDohrSubstruct<TVar>::LoadWeightedResidual(const TPZFMatrix<TVar> &r_global){
80  int i;
81  fLocalWeightedResidual.Resize(fNEquations,1);
82  int neqs = fGlobalEqs.NElements();
83  for (i=0;i<neqs;i++)
84  {
85  std::pair<int,int> ind = fGlobalEqs[i];
86  fLocalWeightedResidual(ind.first,0) = fWeights[ind.first]*r_global.GetVal(ind.second,0);
87  }
88 }
89 
90 template<class TVar>
91 void TPZDohrSubstruct<TVar>::Contribute_rc(TPZFMatrix<TVar> &rc) {
92  int i;
93  TPZFMatrix<TVar> temp;
94  fPhiC.Multiply(fLocalWeightedResidual, temp, 1);
95  for (i=0;i<fCoarseIndex.NElements();i++) {
96  rc(fCoarseIndex[i],0) += temp(i,0);
97  }
98 }
99 
104 template<class TVar>
105 void TPZDohrSubstruct<TVar>::Contribute_rc_local(TPZFMatrix<TVar> &residual_local, TPZFMatrix<TVar> &rc_local) const
106 {
107  fPhiC_Weighted_Condensed.Multiply(residual_local, rc_local, 1);
108 }
109 
110 
111 template<class TVar>
112 void TPZDohrSubstruct<TVar>::Contribute_Kc(TPZMatrix<TVar> &Kc, TPZVec<int> &coaseindex) {
113  int i;
114  int j;
115  for (i=0;i<fCoarseIndex.NElements();i++) {
116  for (j=0;j<fCoarseIndex.NElements();j++) {
117  Kc(fCoarseIndex[i],fCoarseIndex[j]) += fKCi(i,j);
118  }
119  }
120 }
121 
122 template<class TVar>
123 void TPZDohrSubstruct<TVar>::Contribute_v1(TPZFMatrix<TVar> &v1, TPZFMatrix<TVar> &invKc_rc) {
124  int i;
125  // int j;
126  TPZFMatrix<TVar> temp(fCoarseIndex.NElements(), 1);
127  TPZFMatrix<TVar> temp2;
128  //temp = R(ci)*K(c)_inverted*r(c)
129  for (i=0;i<fCoarseIndex.NElements();i++) {
130  temp(i, 0) = invKc_rc(fCoarseIndex[i],0);
131  }
132  //temp2 = Phi*temp
133  fPhiC.Multiply(temp, temp2, 0);
134  //v1 += R(i)_transposta*W(i)*temp2
135 #ifdef ZERO_INTERNAL_RESIDU
136  for(i=0; i<fInternalEqs.NElements(); i++)
137  {
138  temp2(fInternalEqs[i],0) = 0.;
139  }
140 #endif
141  int neqs = fGlobalEqs.NElements();
142  for (i=0;i<neqs;i++) {
143  std::pair<int,int> ind = fGlobalEqs[i];
144  v1(ind.second,0) += fWeights[ind.first]*temp2(ind.first,0);
145  }
146 }
147 
148 template<class TVar>
149 void TPZDohrSubstruct<TVar>::Contribute_v1_local(TPZFMatrix<TVar> &v1_local, TPZFMatrix<TVar> &invKc_rc_local) const {
150  int neqs = fGlobalEqs.NElements();
151  v1_local.Resize(neqs, 1);
152  fPhiC_Weighted_Condensed.Multiply(invKc_rc_local,v1_local);
153 }
154 
155 template<class TVar>
156 void TPZDohrSubstruct<TVar>::Contribute_v2(TPZFMatrix<TVar> &v2) {
157  int i;
158  SolveSystemZi();
159 #ifdef ZERO_INTERNAL_RESIDU
160  for(i=0; i<fInternalEqs.NElements(); i++)
161  {
162  fzi(fInternalEqs[i],0) = 0.;
163  }
164 #endif
165  int neqs = fGlobalEqs.NElements();
166  for (i=0;i<neqs;i++)
167  {
168  std::pair<int,int> ind = fGlobalEqs[i];
169  v2(ind.second,0) += fWeights[ind.first] * fzi(ind.first,0);
170  }
171 }
172 
173 /*
174  * It computes the local contribution to v2.
175  */
176 template<class TVar>
177 void TPZDohrSubstruct<TVar>::Contribute_v2_local(TPZFMatrix<TVar> &residual_local, TPZFMatrix<TVar> &v2_local)
178 {
179  int ncols = residual_local.Cols();
180  TPZFMatrix<TVar> LocalWeightedResidual(fNEquations,ncols,0.);
181  int neqs = fGlobalEqs.NElements();
182  int i;
183  for (int ic=0; ic<ncols; ic++)
184  {
185  for (i=0;i<neqs;i++)
186  {
187  std::pair<int,int> ind = fGlobalEqs[i];
188  LocalWeightedResidual(ind.first,ic) += fWeights[ind.first] * residual_local(i,ic);
189  }
190  }
191  int ncoarse = fCoarseIndex.NElements();
192  // size of the kernel
193  int nnull = fNullPivots.Rows();
194  // number of global indices
195  int nglob = fNEquations;
197  //Constructing I star is the same I star for Phi
198  //C star is the same C star for Phi
199  //Constructing I_lambda
200  TPZFMatrix<TVar> I_lambda(ncoarse+nnull,ncols);
201  I_lambda.Zero();
202  //K_star_inv*C_star_trans is the same used for Phi
203  /* Computing K_star_inv*W(i)*R(i)*r = K_star_inv*fLocalWeightedResidual */
204  TPZFMatrix<TVar> KWeightedResidual(nglob,ncols);
205  fInvertedStiffness.Solve(LocalWeightedResidual,KWeightedResidual);
206  //Obtaining lambda_star
207  TPZFMatrix<TVar> Lambda_star(ncoarse+nnull,1);
208  TPZFMatrix<TVar> CstarKW(ncoarse+nnull,ncols);
209  fC_star.MultAdd(KWeightedResidual,KWeightedResidual,CstarKW,-1,0,0);
210  TPZFMatrix<TVar> temp2(ncoarse+nnull,ncols);
211  I_lambda.Add(CstarKW,temp2);
212  finv.Solve(temp2, Lambda_star);
213 #ifdef LOG4CXX
214  if(logger->isDebugEnabled())
215  {
216  std::stringstream sout;
217  Lambda_star.Print("Lambda_star ",sout);
218  LOGPZ_DEBUG(logger,sout.str())
219  }
220 #endif
221  //Obtaining z(i)
222  TPZFMatrix<TVar> zi(nglob,ncols);
223  temp2.Resize(nglob,ncols);
224  fKeC_star.Multiply(Lambda_star,temp2);
225  temp2 *= -1.;
226  temp2.Add(KWeightedResidual,zi);
227 #ifdef LOG4CXX
228  if(logger->isDebugEnabled())
229  {
230  std::stringstream sout;
231  zi.Print("zi ",sout);
232  LOGPZ_DEBUG(logger,sout.str())
233  }
234 #endif
235 #ifdef ZERO_INTERNAL_RESIDU
236  for(i=0; i<fInternalEqs.NElements(); i++)
237  {
238  zi(fInternalEqs[i],0) = 0.;
239  }
240 #endif
241  v2_local.Resize(neqs, ncols);
242  for (int ic=0; ic<ncols; ic++)
243  {
244  for (i=0;i<neqs;i++)
245  {
246  std::pair<int,int> ind = fGlobalEqs[i];
247  v2_local(i,ic) = fWeights[ind.first] * zi(ind.first,ic);
248  }
249  }
250 
251 
252 }
253 
254 template<class TVar>
255 void TPZDohrSubstruct<TVar>::Contribute_v3(TPZFMatrix<TVar> &v3, const TPZFMatrix<TVar> &r, TPZFMatrix<TVar> &v1Plusv2) const {
256  TPZFNMatrix<100,TVar> vec_t(fNEquations,1,0.);
257  TPZFNMatrix<100,TVar> vec_t2(fNEquations,1,0.);
258  TPZFNMatrix<100,TVar> vec_t3(fInternalEqs.NElements(),1,0.);
259  TPZFNMatrix<100,TVar> inv_sys(fInternalEqs.NElements(),1,0.);
260  int i;
261  int neqs = fGlobalEqs.NElements();
262  //vec_t=R(i)(v1+v2)
263  for (i=0;i<neqs;i++) {
264  std::pair<int,int> ind = fGlobalEqs[i];
265  vec_t(ind.first,0) = v1Plusv2(ind.second,0);
266  }
267  //vec_t2=K(i)*vec_t
268  fStiffness->Multiply(vec_t, vec_t2, 0);
269  //vec_t=R(i)r
270  for (i=0;i<neqs;i++) {
271  std::pair<int,int> ind = fGlobalEqs[i];
272  vec_t(ind.first,0) = r.GetVal(ind.second,0);
273  }
274  //vec_t3=R(Ii)*(vec_t - vec_t2)
275  for (i=0;i<fInternalEqs.NElements();i++) {
276  vec_t3(i,0) = vec_t(fInternalEqs[i],0) - vec_t2(fInternalEqs[i],0);
277  }
278  //inv_sys = temp_inverted * vec_t3
279  fInvertedInternalStiffness.Solve(vec_t3, inv_sys);
280  //vec_t=R(Ii)_transposed * inv_sys
281  vec_t.Zero();
282  for (i=0;i<fInternalEqs.NElements();i++) {
283  vec_t(fInternalEqs[i],0) = inv_sys(i,0);
284  }
285 #ifndef MAKEINTERNAL
286  //v3=v3+R(i)_transposed*vec_t
287  for (i=0;i<neqs;i++) {
288  std::pair<int,int> ind = fGlobalEqs[i];
289  v3(ind.second,0) += vec_t(ind.first,0);
290  }
291 #else
292  v3 = vec_t;
293 #endif
294 
295 }
296 
297 template<class TVar>
298 void TPZDohrSubstruct<TVar>::Contribute_v3_local(TPZFMatrix<TVar> &v1Plusv2, TPZFMatrix<TVar> &v3) const {
299  TPZFNMatrix<100,TVar> vec_t2(fNEquations,1,0.);
300  TPZFNMatrix<100,TVar> vec_t3(fInternalEqs.NElements(),1,0.);
301  TPZFNMatrix<100,TVar> inv_sys(fInternalEqs.NElements(),1,0.);
302  int i;
303  //vec_t2=K(i)*vec_t
304  fStiffness->Multiply(v1Plusv2, vec_t2);
305  for (i=0;i<fInternalEqs.NElements();i++) {
306  vec_t3(i,0) = vec_t2(fInternalEqs[i],0);
307  }
308  //inv_sys = temp_inverted * vec_t3
309  fInvertedInternalStiffness.Solve(vec_t3, inv_sys);
310  //vec_t=R(Ii)_transposed * inv_sys
311  v3.Redim(fNEquations,1);
312  for (i=0;i<fInternalEqs.NElements();i++) {
313  v3(fInternalEqs[i],0) = -inv_sys(i,0);
314  }
315 }
316 
317 template<class TVar>
318 void TPZDohrSubstruct<TVar>::Print(std::ostream &out) const
319 {
320  out << "++++++++++++++++++++++++++++++++++++" << std::endl;
321  out << "Coarse Index " << fCoarseIndex << std::endl;
322  out << "Global Index " << fGlobalIndex << std::endl;
323  out << "Internal Nodes " << fInternalEqs << std::endl;
324  out << "Coarse Nodes " << fCoarseNodes << std::endl;
325  //fEigenVectors.Print("Eigen vectors",out);
326  fC_star.Print("fC_star",out);
327  fKeC_star.Print("fKeC_star", out);
328  fNullPivots.Print("fNullPivots", out);
329  out << "fNEquations = " << fNEquations << std::endl;
330  out << "fCoarseNodes " << fCoarseNodes << std::endl;
331 
332  out << "fCoarseIndex " << fCoarseIndex << std::endl;
333 
334  out << "fGlobalEqs " << fGlobalEqs << std::endl;
335  out << "fInternalEqs " << fInternalEqs << std::endl;
336  out << "fBoundaryEqs " << fBoundaryEqs << std::endl;
337  fLocalLoad.Print("fLocalLoad",out);
338  fLocalWeightedResidual.Print("fLocalWeightedResidual",out);
339  fzi.Print("fzi", out);
340  fAdjustSolution.Print("fAdjustSolution", out);
341 
342 
343  fKCi.Print("Coarse Matrix",out);
344  fStiffness->Print("Stiffness Matrix",out,EMathematicaInput);
345  fInvertedStiffness.Matrix()->Print("Matrix Ke",out,EMathematicaInput);
346  fC.Print("Matrix Ci data structure fC",out,EMathematicaInput);
347  //fEigenVectors.Print("Eigenvectors",out,EMathematicaInput);
348  fPhiC.Print("fPhiC = ",out,EMathematicaInput);
349  out << "fWeights = " << fWeights << endl;
350  fPhiC_Weighted_Condensed.Print("fPhiC_Weighted_Condensed = ", out);
351 }
352 
353 template<class TVar>
354 void TPZDohrSubstruct<TVar>::SolveSystemPhi() {
355  int ncoarse = fCoarseIndex.NElements();
356  int i;
357  // I_star.Print("fIStar = ",out,EMathematicaInput);
358  // std::cout << "Nci: ";
359  // std::cout << ncoarse;// << endl;
360  // std::cout << endl;
361  //Constructing I_lambda
362  TPZFMatrix<TVar> I_lambda(ncoarse+fNullPivots.Rows(),ncoarse);
363  I_lambda.Zero();
364  for (i=0;i<ncoarse;i++) {
365  I_lambda(i,i)=1;
366  }
367  // I_lambda.Print("ILambda = ",out,EMathematicaInput);
368 
369  //Obtaining lambda_star
370  TPZFMatrix<TVar> Lambda_star(ncoarse+fNullPivots.Rows(),ncoarse);
371  // temp1->Print("temp1 = ",out,EMathematicaInput);
372  // out.flush();
373  finv.Solve(I_lambda, Lambda_star);
374 
375 #ifdef LOG4CXX
376  if(logger->isDebugEnabled())
377  {
378  std::stringstream sout;
379  Lambda_star.Print("matrix lambda star",sout );
380  LOGPZ_DEBUG(logger,sout.str())
381  }
382 #endif
383  //Obtaining Phi
385  /*
386  TPZFMatrix<TVar> temp2(fNEquations,ncoarse);
387  C_star.MultAdd(Lambda_star,Lambda_star,temp2,-1,0,1,1);
388  */
389  fPhiC.Resize(fNEquations,ncoarse);
390  fKeC_star.Multiply(Lambda_star,fPhiC);
391  fPhiC *= -1.;
392  ComputeCoarseStiffness();
393  // fPhiC.Print("PhiC = ",out,EMathematicaInput);
394  // fC.Print("Ci = ",std::cout,EMathematicaInput);
395 }
396 
397 template<class TVar>
398 void TPZDohrSubstruct<TVar>::SolveSystemZi() {
399  int ncoarse = fCoarseIndex.NElements();
400  // size of the kernel
401  int nnull = fNullPivots.Rows();
402  // number of global indices
403  int nglob = fNEquations;
404  /* Solving the system for zi */
405  //Constructing I star is the same I star for Phi
406  //C star is the same C star for Phi
407  //Constructing I_lambda
408  TPZFMatrix<TVar> I_lambda(ncoarse+nnull,1);
409  I_lambda.Zero();
410  //K_star_inv*C_star_trans is the same used for Phi
411  /* Computing K_star_inv*W(i)*R(i)*r = K_star_inv*fLocalWeightedResidual */
412  TPZFMatrix<TVar> KWeightedResidual(nglob,1);
413  fInvertedStiffness.Solve(fLocalWeightedResidual,KWeightedResidual);
414  //Obtaining lambda_star
415  TPZFMatrix<TVar> Lambda_star(ncoarse+nnull,1);
416  TPZFMatrix<TVar> CstarKW(ncoarse+nnull,1);
417  fC_star.MultAdd(KWeightedResidual,KWeightedResidual,CstarKW,-1,0,0);
418  TPZFMatrix<TVar> temp2(ncoarse+nnull,1);
419  I_lambda.Add(CstarKW,temp2);
420  finv.Solve(temp2, Lambda_star);
421  //Obtaining z(i)
422  fzi.Resize(nglob,1);
423  temp2.Resize(nglob,1);
424  fKeC_star.Multiply(Lambda_star,temp2);
425  temp2 *= -1.;
426  temp2.Add(KWeightedResidual,fzi);
427 }
428 
429 template<class TVar>
430 void TPZDohrSubstruct<TVar>::ComputeCoarseStiffness() {
431  TPZFMatrix<TVar> temp1(fNEquations,fCoarseIndex.NElements());
432  fKCi.Resize(fCoarseIndex.NElements(),fCoarseIndex.NElements());
433  fStiffness->Multiply(fPhiC,temp1);
434 #ifdef LOG4CXX
435  if(logger->isDebugEnabled())
436  {
437  std::stringstream sout;
438  temp1.Print("stiffness matrix times phi",sout);
439  LOGPZ_DEBUG(logger,sout.str())
440  }
441 #endif
442  fPhiC.MultAdd(temp1,temp1,fKCi,1,0,1);
443 #ifdef LOG4CXX
444  if(logger->isDebugEnabled())
445  {
446  std::stringstream sout;
447  fKCi.Print("Coarse stiffness matrix",sout);
448  LOGPZ_DEBUG(logger,sout.str())
449  }
450 #endif
451 }
452 
453 template<class TVar>
454 void TPZDohrSubstruct<TVar>::ContributeGlobalDiagonal(TPZFMatrix<TVar> &StiffnessDiag) {
455  int i;
456  fWeights.Resize(fNEquations);
457  //fWeights = diag(kci) if u(i) is on coarse or diag(Ki)
458  for (i=0;i<fNEquations;i++) {
459  fWeights[i] = fStiffness->operator()(i,i);
460  }
461  if(fWeightType == CorrectWeight)
462  {
463  for (i=0;i<fNEquations;i++) {
464  fWeights[i] = fStiffness->operator()(i,i);
465  }
466  for (i=0;i<fCoarseNodes.NElements();i++) {
467  fWeights[fCoarseNodes[i]] = fKCi(i,i);
468  }
469  } else {
470  for (i=0;i<fNEquations;i++) {
471  fWeights[i] = 1.;
472  }
473  for (i=0;i<fCoarseNodes.NElements();i++) {
474  fWeights[fCoarseNodes[i]] = 1.;
475  }
476  }
477 #ifdef LOG4CXX
478  {
479  std::stringstream sout;
480  sout << "Weight used for assembly" << fWeights;
481  LOGPZ_DEBUG(logger,sout.str())
482  }
483 #endif
484  //
485  int neqs = fGlobalEqs.NElements();
486  for (i=0;i<neqs;i++) {
487  std::pair<int,int> ind = fGlobalEqs[i];
488  StiffnessDiag(ind.second,0) += fWeights[ind.first];
489  }
490 }
491 
492 template<class TVar>
493 void TPZDohrSubstruct<TVar>::ContributeDiagonalLocal(TPZFMatrix<TVar> &StiffnessDiagLocal) {
494  int i;
495  fWeights.Resize(fNEquations);
496  //fWeights = diag(kci) if u(i) is on coarse or diag(Ki)
497  for (i=0;i<fNEquations;i++) {
498  fWeights[i] = fStiffness->operator()(i,i);
499  }
500  if(fWeightType == CorrectWeight)
501  {
502  for (i=0;i<fNEquations;i++) {
503  fWeights[i] = fStiffness->operator()(i,i);
504  }
505  for (i=0;i<fCoarseNodes.NElements();i++) {
506  fWeights[fCoarseNodes[i]] = fKCi(i,i);
507  }
508  } else {
509  for (i=0;i<fNEquations;i++) {
510  fWeights[i] = 1.;
511  }
512  for (i=0;i<fCoarseNodes.NElements();i++) {
513  fWeights[fCoarseNodes[i]] = 1.;
514  }
515  }
516 #ifdef LOG4CXX
517  {
518  std::stringstream sout;
519  sout << "Weight used for assembly" << fWeights;
520  LOGPZ_DEBUG(logger,sout.str())
521  }
522 #endif
523 
524  int neqs = fGlobalEqs.NElements();
525  StiffnessDiagLocal.Resize(neqs,1);
526  for (i=0;i<neqs;i++) {
527  std::pair<int,int> ind = fGlobalEqs[i];
528  StiffnessDiagLocal(i,0) = fWeights[ind.first];
529  }
530 }
531 
532 template<class TVar>
533 void TPZDohrSubstruct<TVar>::ComputeWeights(TPZFMatrix<TVar> &StiffnessDiag) {
534  int i;
535  //fWeights.Fill(1.);
536  int neqs = fGlobalEqs.NElements();
537  for (i=0;i<neqs;i++) {
538  std::pair<int,int> ind = fGlobalEqs[i];
539  fWeights[ind.first] = fWeights[ind.first] / StiffnessDiag(ind.second,0);
540  }
541  for(i=0; i<fInternalEqs.NElements(); i++)
542  {
543  fWeights[fInternalEqs[i]] = 1.;
544  }
545 #ifdef LOG4CXX
546  {
547  std::stringstream sout;
548  sout << "Weights = " << fWeights;
549  LOGPZ_DEBUG(logger,sout.str())
550  }
551 #endif
552 }
553 
554 template<class TVar>
555 void TPZDohrSubstruct<TVar>::ComputeWeightsLocal(TPZFMatrix<TVar> &StiffnessDiagLocal) {
556  int i;
557  int neqs = fGlobalEqs.NElements();
558  for (i=0;i<neqs;i++) {
559  std::pair<int,int> ind = fGlobalEqs[i];
560  fWeights[ind.first] = fWeights[ind.first] / StiffnessDiagLocal(i,0);
561  }
562  for(i=0; i<fInternalEqs.NElements(); i++)
563  {
564  fWeights[fInternalEqs[i]] = 1.;
565  }
566 #ifdef LOG4CXX
567  {
568  std::stringstream sout;
569  sout << "Weights = " << fWeights;
570  LOGPZ_DEBUG(logger,sout.str())
571  }
572 #endif
573  int c,nc = fPhiC.Cols();
574 #ifdef ZERO_INTERNAL_RESIDU
575  // eliminate the values of v2 for internal equations
576  int nint = fInternalEqs.NElements();
577  for(i=0; i<nint; i++)
578  {
579  for(c=0; c<nc; c++)
580  {
581  fPhiC(fInternalEqs[i],c) = 0.;
582  }
583  }
584 #endif
585  fPhiC_Weighted_Condensed.Resize(neqs,fPhiC.Cols());
586  for (i=0;i<neqs;i++)
587  {
588  for(c=0; c<nc; c++)
589  {
590  std::pair<int,int> ind = fGlobalEqs[i];
591  fPhiC_Weighted_Condensed(i,c) = fPhiC(ind.first,c)*fWeights[ind.first];
592  }
593  }
594 }
595 
596 template<class TVar>
597 void TPZDohrSubstruct<TVar>::ContributeKU(const TVar alpha, const TPZFMatrix<TVar> &uglobal, TPZFMatrix<TVar> &z) const
598 {
599  int i,j;
600  int nglob = fNEquations;
601  int nglobglob = uglobal.Rows();
602  int ncols = uglobal.Cols();
603  TPZFMatrix<TVar> uloc(nglob,ncols);
604  TPZFNMatrix<100,TVar> temp1(nglob,ncols),v3(nglob,ncols);
605  TPZFNMatrix<100,TVar> temp1glob(nglobglob,ncols,0.),zero(nglobglob,ncols,0.),v3glob(nglobglob,ncols,0.);
606  TPZFMatrix<TVar> kuloc;
607  uloc.Zero();
608  int neqs = fGlobalEqs.NElements();
609  /* Performing R(i)*u=uloc */
610  for (i=0;i<neqs;i++) {
611  std::pair<int,int> ind = fGlobalEqs[i];
612  for(j=0; j<ncols; j++)
613  {
614  uloc(ind.first,j) = uglobal.Get(ind.second,j);
615  }
616  }
617 #ifdef LOG4CXX
618  {
619  std::stringstream sout;
620  sout << "Value of the local solution = ";
621  uloc.Print("uloc " ,sout);
622  LOGPZ_DEBUG(logger,sout.str())
623  }
624 #endif
625 #ifdef ZERO_INTERNAL_RESIDU
626  // convert the local temp1 matrix para global
627  for (i=0;i<neqs;i++) {
628  for(j=0; j<ncols; j++)
629  {
630  std::pair<int,int> ind = fGlobalEqs[i];
631  temp1glob(ind.second,j) = uloc(ind.first,j);
632  }
633  }
634  // zero out the internal residu
635  this->Contribute_v3(v3glob,zero,temp1glob);
636 #ifndef MAKEINTERNAL
637  // convert the local temp1 matrix para global
638  for (i=0;i<neqs;i++) {
639  for(j=0; j<ncols; j++)
640  {
641  std::pair<int,int> ind = fGlobalEqs[i];
642  v3(ind.first,j) = v3glob(ind.second,j);
643  }
644  }
645 #else
646  v3 = v3glob;
647 #endif
648 #ifdef LOG4CXX
649  {
650  std::stringstream sout;
651  sout << "Contribution of v3";
652  v3.Print("v3",sout);
653  LOGPZ_DEBUG(logger,sout.str())
654  }
655 #endif
656  uloc += v3;
657 #ifdef LOG4CXX
658  {
659  std::stringstream sout;
660  sout << "New value of uloc";
661  uloc.Print("uloc",sout);
662  LOGPZ_DEBUG(logger,sout.str())
663  }
664 #endif
665 #endif
666  /* Computing K(i)*uloc and storing in temp1 */
667  temp1.Resize(nglob,ncols);
668  fStiffness->Multiply(uloc,temp1);
669 #ifdef LOG4CXX
670  {
671  std::stringstream sout;
672  sout << "After multiplying the solution temp1 = ";
673  temp1.Print("temp1 " ,sout);
674  LOGPZ_DEBUG(logger,sout.str())
675  }
676 #endif
677  int zcols = z.Cols();
678  for (i=0;i<neqs;i++) {
679  std::pair<int,int> ind = fGlobalEqs[i];
680  /* Sum row "i" of temp1 with row "fGlobalNodes[i]" of z */
681  for (j=0;j<zcols;j++) {
682  z(ind.second,j) += alpha*temp1(ind.first,j);
683  }
684  }
685 }
686 
687 template<class TVar>
688 void TPZDohrSubstruct<TVar>::ContributeKULocal(const TVar alpha, const TPZFMatrix<TVar> &u, TPZFMatrix<TVar> &z) const
689 {
690  int i,j;
691  int nglob = fNEquations;
692  int ncols = u.Cols();
693  TPZFMatrix<TVar> uloc(nglob,ncols,0.);
694  TPZFNMatrix<100,TVar> temp1(nglob,ncols),v3(nglob,ncols);
695  TPZFMatrix<TVar> kuloc;
696  uloc.Zero();
697  int neqs = u.Rows();
698  /* Performing R(i)*u=uloc */
699  for (i=0;i<neqs;i++) {
700  std::pair<int,int> ind = fGlobalEqs[i];
701  for(j=0; j<ncols; j++)
702  {
703  uloc(ind.first,j) = u.Get(i,j);
704  }
705  }
706 #ifdef LOG4CXX
707  {
708  std::stringstream sout;
709  sout << "Value of the local solution = ";
710  uloc.Print("uloc " ,sout);
711  LOGPZ_DEBUG(logger,sout.str())
712  }
713 #endif
714  // *****************************************************************
715 #ifdef ZERO_INTERNAL_RESIDU
716  // zero out the internal residu
717  this->Contribute_v3_local(uloc,v3);
718 #ifdef LOG4CXX
719  {
720  std::stringstream sout;
721  sout << "Contribution of v3";
722  v3.Print("v3",sout);
723  LOGPZ_DEBUG(logger,sout.str())
724  }
725 #endif
726  uloc += v3;
727 #ifdef LOG4CXX
728  {
729  std::stringstream sout;
730  sout << "New value of uloc";
731  uloc.Print("uloc",sout);
732  LOGPZ_DEBUG(logger,sout.str())
733  }
734 #endif
735 #endif
736  /* Computing K(i)*uloc and storing in temp1 */
737  temp1.Resize(nglob,ncols);
738  fStiffness->Multiply(uloc,temp1);
739 #ifdef LOG4CXX
740  {
741  std::stringstream sout;
742  sout << "After multiplying the solution temp1 = ";
743  temp1.Print("temp1 " ,sout);
744  LOGPZ_DEBUG(logger,sout.str())
745  }
746 #endif
747  int zcols = z.Cols();
748  for (i=0;i<neqs;i++) {
749  std::pair<int,int> ind = fGlobalEqs[i];
750  /* Sum row "i" of temp1 with row "fGlobalNodes[i]" of z */
751  for (j=0;j<zcols;j++) {
752  z(i,j) += alpha*temp1(ind.first,j);
753  }
754  }
755 }
756 
757 template<class TVar>
758 void TPZDohrSubstruct<TVar>::Initialize() {
759  PrepareSystems();
760  SolveSystemPhi();
761  ComputeCoarseStiffness();
762 
763 }
764 
765 template<class TVar>
766 void TPZDohrSubstruct<TVar>::PrepareSystems() {
767  int ncoarse = fCoarseIndex.NElements();
768  //TPZFMatrix<TVar> C_transposed;
769  /****************************************
770  * TEMPORARIO!!!!!!!!!!!!!!
771  */
772  // int nc = fC.Cols();
773  // fC.Redim(1,nc);
774  // fC(0,nc-1) = 1.;
775 
776  /*
777  *******************************************/
780  // ofstream out("saida.nb");
781  // this->fStiffness->Print("fK = ",out,EMathematicaInput);
782  // fC.Print("fCi = ",out,EMathematicaInput);
783  TPZFMatrix<TVar> C_transposed(fC.Cols(),fC.Rows());
784  fC.Transpose(&C_transposed);
785  //Ke_inv*C(i)_trans
786  TPZFMatrix<TVar> KeC(fNEquations,ncoarse);
787  fInvertedStiffness.Solve(C_transposed,KeC);
789 #ifdef LOG4CXX
790  {
791  std::stringstream sout;
792  sout << "Number of zero pivots of the substructure : " << fInvertedStiffness.Singular().size();
793  LOGPZ_DEBUG(logger,sout.str());
794  }
795 #endif
796  // std::cout << KeC;
797  // int neq = fInvertedStiffness.Matrix()->Rows();
798  //Constructing C barra = NullPivots
799  fNullPivots.Resize(fInvertedStiffness.Singular().size(),fInvertedStiffness.Matrix()->Cols());
800  fNullPivots.Zero();
801  // fEigenVectors.Redim(fInvertedStiffness.Matrix()->Rows(),fInvertedStiffness.Singular().size());
802  if(fNullPivots.Rows())
803  {
804  // out << "NullPivots = {";
805  int count = 0;
806  std::list<int64_t>::iterator it = fInvertedStiffness.Singular().begin();
807  for(;it != fInvertedStiffness.Singular().end(); it++,count++)
808  {
809  // out << *it;
810  // if(count < fInvertedStiffness.Singular().size()-1) out << ",";
811  fNullPivots(count,*it) = 1.;
812  }
813  // out << "};\n";
814  // fInvertedStiffness.Solve(fEigenVectors,fEigenVectors);
815  }
816  // NullPivots.Print("fNullPivots = ",out,EMathematicaInput);
817 
818  //Constructing I star
819  TPZFMatrix<TVar> I_star(ncoarse+fNullPivots.Rows(),ncoarse+fNullPivots.Rows());
820  I_star.Zero();
821  int i;
822  for (i=ncoarse;i<ncoarse+fNullPivots.Rows();i++) {
823  I_star(i,i)=1;
824  }
825  //Constructing C_star and storing in fC_star
826  fC_star.Resize(ncoarse+fNullPivots.Rows(),fNEquations);
827  fC_star.PutSub(0,0,fC);
828  fC_star.PutSub(ncoarse,0,fNullPivots);
829  // C_star.Print("CStar = ",out,EMathematicaInput);
830  //constructing K_star_inv*C_star_trans and storing it in fKeC_star
831  fKeC_star.Resize(fNEquations,ncoarse+fNullPivots.Rows());
832  TPZFMatrix<TVar> C_star_trans(fNEquations,ncoarse+fNullPivots.Rows());
833  fC_star.Transpose(&C_star_trans);
834  fInvertedStiffness.Solve(C_star_trans,fKeC_star);
835  //Constructing StepSolver finv
836  TPZFMatrix<TVar> *temp1 = new TPZFMatrix<TVar>(ncoarse+fNullPivots.Rows(),ncoarse+fNullPivots.Rows());
837  fC_star.MultAdd(fKeC_star,I_star,*temp1,-1,1,0);
838  finv.SetMatrix(temp1);
839  finv.SetDirect(ELU);
840 }
841 
846 template<class TVar>
847 void TPZDohrSubstruct<TVar>::AdjustResidual(TPZFMatrix<TVar> &r_global)
848 {
849  int nglob = fNEquations;
850  int nint = this->fInternalEqs.NElements();
851  TPZFMatrix<TVar> rint(nint,1,0.),radapt(nglob,1,0.);
852  int i;
853  for(i=0; i<nint; i++)
854  {
855  if(fGlobalIndex[fInternalEqs[i]] >= 0)
856  {
857  rint(i,0) = r_global(fGlobalIndex[fInternalEqs[i]]);
858  }
859  }
860 
861  // initializar radapt com zero
862  int neqs = fGlobalEqs.NElements();
863  for(i=0; i<neqs; i++)
864  {
865  std::pair<int,int> ind = fGlobalEqs[i];
866  radapt(ind.first,0) = r_global(ind.second,0);
867  }
868  fInvertedInternalStiffness.Solve(rint,fAdjustSolution);
869 #ifdef LOG4CXX
870  {
871  std::stringstream sout;
872  sout << "Adjusted internal solution";
873  fAdjustSolution.Print("fAdjustSolution",sout);
874  rint.Print("Internal residual",sout);
875  LOGPZ_DEBUG(logger,sout.str())
876  }
877 #endif
878  TPZFMatrix<TVar> rloc(nglob,1,0.),radjust(nglob,1,0.);
879  for(i=0; i<nint; i++)
880  {
881  rloc(fInternalEqs[i],0) = fAdjustSolution(i,0);
882  }
883  fStiffness->MultAdd(rloc,radapt,radjust,-1.,1.);
884  for(i=0; i<nint; i++)
885  {
886  if(fabs(radjust(fInternalEqs[i],0)) >= 1.e-10)
887  {
888  std::cout << "Internal node " << i << " was not zeroed " << radjust(fInternalEqs[i],0) << std::endl;
889  }
890  radjust(fInternalEqs[i],0) = 0.;
891  }
892 
893  // somar a contribuicao no r_global (tarefa)
894  for(i=0; i<neqs; i++)
895  {
896  std::pair<int,int> ind = fGlobalEqs[i];
897  r_global(ind.second,0) = radjust(ind.first,0);
898  }
899 }
900 
904 template<class TVar>
905 void TPZDohrSubstruct<TVar>::AddInternalSolution(TPZFMatrix<TVar> &sol)
906 {
907  int nglob = fNEquations;
908  int nint = this->fInternalEqs.NElements();
909  TPZFMatrix<TVar> rint(nint,1),radapt(nglob,1,0.);
910  int i,c,nc = sol.Cols();
911  for(i=0; i<nint; i++) for(c=0; c<nc; c++)
912  {
913  if(fGlobalIndex[fInternalEqs[i]] >=0) sol(fGlobalIndex[fInternalEqs[i]],c) += this->fAdjustSolution(i,c);
914  }
915 
916 }
917 
918 template class TPZDohrSubstruct<float>;
919 template class TPZDohrSubstruct<double>;
920 template class TPZDohrSubstruct<long double>;
921 
922 //#ifdef STATE_COMPLEX
923 template class TPZDohrSubstruct<std::complex<float> >;
924 template class TPZDohrSubstruct<std::complex<double> >;
925 template class TPZDohrSubstruct<std::complex<long double> >;
926 //#endif
fabs
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ fabs
Definition: tfadfunc.h:140
TPZDohrSubstruct::Contribute_Kc
void Contribute_Kc(TPZMatrix< TVar > &Kc, TPZVec< int > &coarseindex)
It computes the local contribution to K(c)
Definition: tpzdohrsubstruct.cpp:112
TPZDohrSubstruct::Contribute_v1
void Contribute_v1(TPZFMatrix< TVar > &v1, TPZFMatrix< TVar > &invKc_rc)
It computes the local contribution to v1.
Definition: tpzdohrsubstruct.cpp:123
TPZMatrix::PutSub
virtual int PutSub(const int64_t sRow, const int64_t sCol, const TPZFMatrix< TVar > &Source)
It puts submatrix Source on actual matrix structure.
Definition: pzmatrix.cpp:574
pzlog.h
Contains definitions to LOGPZ_DEBUG, LOGPZ_INFO, LOGPZ_WARN, LOGPZ_ERROR and LOGPZ_FATAL, and the implementation of the inline InitializePZLOG(string) function using log4cxx library or not. It must to be called out of "#ifdef LOG4CXX" scope.
TPZDohrSubstruct::Contribute_v3
void Contribute_v3(TPZFMatrix< TVar > &v3, const TPZFMatrix< TVar > &r, TPZFMatrix< TVar > &v1Plusv2) const
It computes the local contribution to v(3)
Definition: tpzdohrsubstruct.cpp:255
TPZDohrSubstruct::ContributeKU
void ContributeKU(const TVar alpha, const TPZFMatrix< TVar > &uglobal, TPZFMatrix< TVar > &z) const
Contribute to the global matrix vector multiplication. It&#39;s needed to the MultAdd method of TPZDohrMa...
Definition: tpzdohrsubstruct.cpp:597
TPZDohrSubstruct::SolveSystemPhi
void SolveSystemPhi()
Solves the system for Phi and for v2.
Definition: tpzdohrsubstruct.cpp:354
TPZDohrSubstruct::ContributeTestV1
void ContributeTestV1(TPZFMatrix< TVar > &testV1, int NumCoarse)
Definition: tpzdohrsubstruct.cpp:33
tpzdohrsubstruct.h
Contains the TPZDohrSubstruct class which implements sub structure matrices using Dohrman algorithm...
TPZDohrSubstruct::PrepareSystems
void PrepareSystems()
It prepares the datas for solving systems for phi and zi.
Definition: tpzdohrsubstruct.cpp:766
TPZDohrSubstruct::ComputeWeights
void ComputeWeights(TPZFMatrix< TVar > &StiffnessDiag)
Computes the weight matrix.
Definition: tpzdohrsubstruct.cpp:533
TPZMatrix::Get
virtual const TVar & Get(const int64_t row, const int64_t col) const
Get value with bound checking.
Definition: pzmatrix.h:855
TPZDohrSubstruct::ContributeResidual
void ContributeResidual(TPZFMatrix< TVar > &u, TPZFMatrix< TVar > &r)
ContributeResidual.
Definition: tpzdohrsubstruct.cpp:61
TPZDohrSubstruct::Contribute_v2
void Contribute_v2(TPZFMatrix< TVar > &v2)
It computes the local contribution to v2.
Definition: tpzdohrsubstruct.cpp:156
TPZDohrSubstruct
Implements sub structure matrices using Dohrman algorithm. Sub Structure.
Definition: tpzdohrsubstruct.h:28
TPZMatrix::GetSub
virtual int GetSub(const int64_t sRow, const int64_t sCol, const int64_t rowSize, const int64_t colSize, TPZFMatrix< TVar > &Target) const
Gets submatrix storing it on Target.
Definition: pzmatrix.cpp:602
TPZDohrSubstruct::ContributeKULocal
void ContributeKULocal(const TVar alpha, const TPZFMatrix< TVar > &u, TPZFMatrix< TVar > &z) const
Compute the multiplication of the local stiffness matrix with the vector u.
Definition: tpzdohrsubstruct.cpp:688
TPZDohrSubstruct::AdjustResidual
void AdjustResidual(TPZFMatrix< TVar > &r_global)
Adjust the residual to reflect a static condensation.
Definition: tpzdohrsubstruct.cpp:847
TPZDohrSubstruct::SolveSystemZi
void SolveSystemZi()
Solves the system for zi.
Definition: tpzdohrsubstruct.cpp:398
TPZFMatrix::Zero
int Zero() override
Makes Zero all the elements.
Definition: pzfmatrix.h:651
TPZDohrSubstruct::Contribute_rc
void Contribute_rc(TPZFMatrix< TVar > &rc)
It computes the local contribution to r(c).
Definition: tpzdohrsubstruct.cpp:91
TPZDohrSubstruct::ComputeCoarseStiffness
void ComputeCoarseStiffness()
Computes K(ci) and stores it in fKCi.
Definition: tpzdohrsubstruct.cpp:430
LOGPZ_DEBUG
#define LOGPZ_DEBUG(A, B)
Define log for debug info.
Definition: pzlog.h:87
TPZMatrix::Add
virtual void Add(const TPZMatrix< TVar > &A, TPZMatrix< TVar > &res) const
It adds itself to TPZMatrix<TVar>A putting the result in res.
Definition: pzmatrix.cpp:64
ELU
Definition: pzmatrix.h:52
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
TPZFMatrix
Full matrix class. Matrix.
Definition: pzfmatrix.h:32
TPZFMatrix::Redim
int Redim(const int64_t newRows, const int64_t newCols) override
Redimension a matrix and ZERO your elements.
Definition: pzfmatrix.h:616
TPZDohrSubstruct::LoadWeightedResidual
void LoadWeightedResidual(const TPZFMatrix< TVar > &r_global)
It computes W(i)*R(i)*r and stores it in fLocalWeightedResidual;.
Definition: tpzdohrsubstruct.cpp:79
TPZDohrSubstruct::Print
void Print(std::ostream &out) const
Definition: tpzdohrsubstruct.cpp:318
TPZDohrSubstruct::Contribute_v1_local
void Contribute_v1_local(TPZFMatrix< TVar > &v1_local, TPZFMatrix< TVar > &invKc_rc) const
It computes the local contribution to v1.
Definition: tpzdohrsubstruct.cpp:149
TPZDohrSubstruct::AddInternalSolution
void AddInternalSolution(TPZFMatrix< TVar > &sol)
Add the internal solution to the final result.
Definition: tpzdohrsubstruct.cpp:905
TPZMatrix::Multiply
virtual void Multiply(const TPZFMatrix< TVar > &A, TPZFMatrix< TVar > &res, int opt=0) const
It mutiplies itself by TPZMatrix<TVar>A putting the result in res.
Definition: pzmatrix.cpp:1916
TPZDohrSubstruct::ContributeDiagonalLocal
void ContributeDiagonalLocal(TPZFMatrix< TVar > &StiffnessDiag)
Computes the contribution of each substructure node to global Stiffness diagonal (or something like t...
Definition: tpzdohrsubstruct.cpp:493
EMathematicaInput
Definition: pzmatrix.h:55
TPZFMatrix::MultAdd
virtual void MultAdd(const TPZFMatrix< TVar > &x, const TPZFMatrix< TVar > &y, TPZFMatrix< TVar > &z, const TVar alpha=1., const TVar beta=0., const int opt=0) const override
It computes z = beta * y + alpha * opt(this)*x but z and x can not overlap in memory.
Definition: pzfmatrix.cpp:757
TPZMatrix::Cols
int64_t Cols() const
Returns number of cols.
Definition: pzmatrix.h:809
TPZMatrix::Print
virtual void Print(std::ostream &out) const
Definition: pzmatrix.h:253
TPZFMatrix::Resize
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Definition: pzfmatrix.cpp:1016
TPZDohrSubstruct::Contribute_v2_local
void Contribute_v2_local(TPZFMatrix< TVar > &residual_local, TPZFMatrix< TVar > &v2_local)
It computes the local contribution to v2.
Definition: tpzdohrsubstruct.cpp:177
TPZDohrSubstruct::Initialize
void Initialize()
Initializes the substructure.
Definition: tpzdohrsubstruct.cpp:758
TPZFMatrix::GetVal
const TVar & GetVal(const int64_t row, const int64_t col) const override
Get values without bounds checking This method is faster than "Get" if DEBUG is defined.
Definition: pzfmatrix.h:566
TPZFMatrix::Transpose
void Transpose(TPZMatrix< TVar > *const T) const override
It makes *T the transpose of current matrix.
Definition: pzfmatrix.cpp:1073
TPZDohrSubstruct::Contribute_v3_local
void Contribute_v3_local(TPZFMatrix< TVar > &v3, TPZFMatrix< TVar > &v1Plusv2) const
It computes the local contribution to v(3)
Definition: tpzdohrsubstruct.cpp:298
TPZFNMatrix
Non abstract class which implements full matrices with preallocated storage with (N+1) entries...
Definition: pzfmatrix.h:716
TPZMatrix
Root matrix class (abstract). Matrix.
Definition: pzmatrix.h:60
TPZDohrSubstruct::Contribute_rc_local
void Contribute_rc_local(TPZFMatrix< TVar > &residual_local, TPZFMatrix< TVar > &rc_local) const
It computes the local contribution to r(c).
Definition: tpzdohrsubstruct.cpp:105
Generated on Sun Aug 16 2020 11:17:11 for NeoPZ by   doxygen 1.8.13