00001 #include <stdio.h>
00002
00003 #include "z2_smoothing.h"
00004 #include "gadaptivity.h"
00005 #include "adaptivity.h"
00006 #include "global.h"
00007 #include "alias.h"
00008 #include "elemhead.h"
00009 #include "vector.h"
00010 #include "matrix.h"
00011 #include "skyline.h"
00012
00013
00014
00015
00016 z2_smoothing::z2_smoothing (long ncomp)
00017 {
00018 z2_smoothing::ncomp = ncomp;
00019 nn = Mt->nn;
00020 ne = Mt->ne;
00021 flags = Ada->give_adaptflag();
00022 }
00023
00024 z2_smoothing::~z2_smoothing () { }
00025
00026
00027
00028 void z2_smoothing::run (vector *rsigfull)
00029 {
00030 fprintf(stdout,"\n ============== Z2_smoothing ===============");
00031
00032 vector ainv,ntdbr;
00033 skyline *ntn_sky;
00034 char label[50];
00035
00036 allocv (ncomp*nn,ntdbr);
00037 compute_ntdbr (ntdbr);
00038
00039 if (1) {
00040 ntn_sky = new skyline;
00041
00042 alloc_ntn (ntn_sky);
00043 compute_ntn_sky (ntn_sky);
00044 compute_rsigfull (ntn_sky,ntdbr,rsigfull);
00045
00046 delete ntn_sky;
00047 }
00048 else {
00049 fprintf (stdout,"\n used diagonalization");
00050 allocv (nn,ainv);
00051 compute_ainv (ainv);
00052 compute_rsigfull (ainv,ntdbr,rsigfull);
00053 }
00054
00055
00056
00057
00058 sprintf(label, "ntdbr:");
00059 fprintf_vector(Out,ntdbr,label,9);
00060
00061
00062
00063
00064 }
00065
00066
00067 void z2_smoothing::compute_ntdbr (vector &ntdbr)
00068 {
00069 long i,j,te,nne;
00070 ivector cn;
00071 vector ntdbri;
00072
00073 for (i=0;i<ne;i++){
00074 te = Mt->give_elem_type (i);
00075 nne = Mt->give_nne (i);
00076 allocv (ncomp*nne,ntdbri);
00077
00078 switch (te){
00079 case planeelementlt:{ Pelt->ntdbr_vector (i,ntdbri); break; }
00080 case planeelementlq:{ Pelq->ntdbr_vector (i,ntdbri); break; }
00081 case planeelementqt:{ Peqt->ntdbr_vector (i,ntdbri); break; }
00082 case planeelementqq:{ Peqq->ntdbr_vector (i,ntdbri); break; }
00083 case lineartet:{ Ltet->ntdbr_vector (i,ntdbri); break; }
00084 default:{
00085 fprintf (stderr,"\n\n unknown element type is required in function");
00086 fprintf (stderr," adaptivity::compute_ntdbr (%s, line %d)",__FILE__,__LINE__);
00087 }
00088 }
00089
00090 allocv (nne,cn);
00091 give_adapt_code_numbers (i,cn);
00092 for (j=0;j<ncomp;j++){
00093 locglob (ntdbr.a + j*nn,ntdbri.a + j*nne,cn.a,nne);
00094 }
00095 destrv (ntdbri); destrv (cn);
00096 }
00097 }
00098
00099 void z2_smoothing::alloc_ntn (skyline *ntn_sky)
00100 {
00101 ntn_sky->allocadr (nn);
00102
00103 column_lengths_nn (ntn_sky);
00104
00105 ntn_sky->addresses ();
00106 ntn_sky->neglobmat ();
00107 ntn_sky->allocglomat ();
00108
00109 if (flags & 1)
00110 fprintf(stdout,"\n number of skyline entries negm %ld",ntn_sky->negm);
00111 }
00112
00113
00114 void z2_smoothing::column_lengths_nn (skyline *ntn_sky)
00115 {
00116 long i,nne;
00117 ivector cn;
00118
00119 for (i=0;i<ne;i++){
00120 nne = Mt->give_nne (i);
00121 allocv (nne,cn);
00122 give_adapt_code_numbers (i,cn);
00123 ntn_sky->column_lengths_elem (cn.a,nne);
00124 destrv (cn);
00125 }
00126 }
00127
00128 void z2_smoothing::give_adapt_code_numbers (long eid,ivector &cn)
00129 {
00130 long i,nne;
00131
00132 nne = cn.n;
00133 Mt->give_elemnodes (eid,cn);
00134
00135 for (i=0;i<nne;i++)
00136 cn[i] += 1;
00137 }
00138
00139 void z2_smoothing::compute_ntn_sky (skyline *ntn_sky)
00140 {
00141 long i,te,nne;
00142 ivector cn;
00143 matrix ntni;
00144
00145 for (i=0;i<ne;i++){
00146 te = Mt->give_elem_type (i);
00147 nne = Mt->give_nne (i);
00148 allocm (nne,nne,ntni);
00149
00150 switch (te){
00151 case planeelementlt:{ Pelt->ntn_matrix (i,ntni); break; }
00152 case planeelementlq:{ Pelq->ntn_matrix (i,ntni); break; }
00153 case planeelementqt:{ Peqt->ntn_matrix (i,ntni); break; }
00154 case planeelementqq:{ Peqq->ntn_matrix (i,ntni); break; }
00155 case lineartet :{ Ltet->ntn_matrix (i,ntni); break; }
00156 default:{
00157 fprintf (stderr,"\n\n unknown element type is required in function");
00158 fprintf (stderr," adaptivity::compute_ntn_sky (%s, line %d)",__FILE__,__LINE__);
00159 }
00160 }
00161
00162 allocv (nne,cn);
00163 give_adapt_code_numbers (i,cn);
00164 ntn_sky->localized (ntni.a,cn.a,nne);
00165 destrm (ntni); destrv(cn);
00166 }
00167 }
00168
00169 void z2_smoothing :: compute_rsigfull (skyline *ntn_sky,vector &ntdbr, vector *rsigfull)
00170 {
00171 long i, j;
00172 vector sol(nn);
00173
00174 ntn_sky->ldl_sky (sol.a,ntdbr.a,Mp->zero,2);
00175
00176 for (i=0;i<ncomp;i++) {
00177 ntn_sky->ldl_sky (sol.a, ntdbr.a + i*nn, Mp->zero, 3);
00178
00179 for (j=0; j<nn; j++)
00180 rsigfull[j].a[i] = sol.a[j];
00181 }
00182 }
00183
00184 void z2_smoothing::compute_ainv (vector &ainv)
00185 {
00186 long i,te,nne;
00187 ivector cn;
00188 vector ntnvi,ntn(nn);
00189 matrix ntnmi;
00190
00191 for (i=0;i<ne;i++){
00192 te = Mt->give_elem_type (i);
00193 nne = Mt->give_nne (i);
00194 allocv (nne,ntnvi);
00195 allocm (nne,nne,ntnmi);
00196
00197 switch (te){
00198 case planeelementlt:{ Pelt->ntn_matrix (i,ntnmi); break; }
00199 case planeelementlq:{ Pelq->ntn_matrix (i,ntnmi); break; }
00200 case planeelementqt:{ Peqt->ntn_matrix (i,ntnmi); break; }
00201 case planeelementqq:{ Peqq->ntn_matrix (i,ntnmi); break; }
00202 case lineartet:{ Ltet->ntn_matrix (i,ntnmi); break; }
00203 default:{
00204 fprintf (stderr,"\n\n unknown element type is required in function");
00205 fprintf (stderr," adaptivity::compute_ainv (%s, line %d)",__FILE__,__LINE__);
00206 }
00207 }
00208
00209 ntnmtov (ntnmi,ntnvi);
00210
00211 allocv (nne,cn);
00212 give_adapt_code_numbers (i,cn);
00213 locglob (ntn.a,ntnvi.a,cn.a,nne);
00214 destrv (ntnvi); destrm (ntnmi); destrv (cn);
00215 }
00216
00217 for (i=0;i<nn;i++)
00218 ainv[i] = 1/ntn[i];
00219 }
00220
00221
00222 void z2_smoothing::compute_rsigfull (vector &ainv,vector &ntdbr, vector *rsigfull)
00223 {
00224 long i,j;
00225
00226 for (i=0; i<ncomp; i++)
00227 for (j=0; j<nn; j++)
00228
00229 rsigfull[j].a[i] = ainv[j] * ntdbr[j + i*nn];
00230 }