Set up of the problem, construction of the main matrices and of the Riesz Representation of right hand side (free from parameters).
Remember \begin{array}{rcl} -\kappa_b\Delta u & =& 0 \quad\mbox{on each block } B_b,\\ -\kappa_b\frac{\partial u}{\partial n}&=& Bi \cdot\,u \quad \mbox{on } \Gamma_b=\partial B_b\cap \partial\Omega,\quad b=1,2,... (Bi>0 \mbox{ Biot number})\\ -\kappa_b\frac{\partial u}{\partial n}&=& q \quad \mbox{on } \Gamma_{root} \end{array}
 |
| Conductivities |
Variational formulation:
Find \(u\in H^1(\Omega)\) such that for all test \(v\in H^1(\Omega)\) \begin{array}{rcl} \kappa_0\int_{C}\nabla u \cdot \nabla v+(\kappa_0+1)\int_{Cleft}\nabla u \cdot \nabla v+(\kappa_0-1)\int_{Cright}\nabla u \cdot \nabla v+ \int_{\partial \Omega\setminus \Gamma_{root}} u \,v& =& \int_{\Gamma_{root}} 10 \,v. \end{array}
Make Riesz Representation of RHS : RieszRepresentation_RHS.idp Info on that..... Riesz
Download here : settingANDproblem.idp Next, go to Offline Computation..... Offline
//
// HERE WE SETUP THE PROBLEM
//
verbosity=0;
load "medit"
//
// In case huge matrices
//load "UMFPACK64"
//
load "msh3"
macro Grad3(u) [dx(u),dy(u),dz(u)] // EOM
//
// read the 3D mesh
//
mesh3 Th3("Th3.mesh");
fespace Vh(Th3,P13d);
cout<<"DOFs "<<Vh.ndof<<endl;
Vh u,v,fh,qh, one=1;
//
// Get volume of computational domain
//
real volu=int3d(Th3)(one);
real invvol=1/volu;
//
// Parameter k0 in [k0min,k0max]
// kleft, kright y Bi funciones de k0
//
real k0,k0min,k0max,kleft,kright,Bi,aa;
include "data.idp"
k0=k0min;
//
// Other parameters depend on k0
//
include "parameters.idp";
//
// Diffusivities
//
fespace Mh(Th3,P03d);
//
Mh kappa0;
kappa0=k0*(region==0)+k0*(region==2)
+k0*(region==4)+k0*(region==6)
+k0*(region==8)+k0*(region==10);
Mh kappa1;
kappa1=kleft*(region==1)+kright*(region==3)
+kleft*(region==5)+kright*(region==7)
+kleft*(region==9)+kright*(region==11);
Mh kappa=kappa0+kappa1;
plot (kappa,wait=0,cmm="KAPPA ",fill=1,nbiso=40,value=1);
//medit("Kappa ",Th3,kappa);
//
// Global problem
//
problem Heatdiffusion(u,v,solver=sparsesolver) =
int3d(Th3,0,2,4,6,8,10)(k0*Grad3(u)' *Grad3(v))
+int3d(Th3,1,5,9)(kleft*Grad3(u)' *Grad3(v))
+int3d(Th3,3,7,11)(kright*Grad3(u)' *Grad3(v))
+int3d(Th3)(aa*u *v) // add elipticidad
+int2d(Th3,10,12,14,16,18,110)(Bi*u*v)// central column front
+int2d(Th3,11,15,19)(Bi*u*v)// front left
+int2d(Th3,13,17,111)(Bi*u*v)// front right
+int2d(Th3,30,32,34,36,38,310)(Bi*u*v)// central column back
+int2d(Th3,31,35,39)(Bi*u*v)// back left
+int2d(Th3,33,37,311)(Bi*u*v)// back right
+int2d(Th3,40,44,48)(Bi*u*v)// left viewpoint central column
+int2d(Th3,41,45,49)(Bi*u*v)// left viewpoint left blocks
+int2d(Th3,20,24,28)(Bi*u*v)// right viewpoint central column
+int2d(Th3,23,27,211)(Bi*u*v)// right viewpoint right blocks
+int2d(Th3,199,201)(Bi*u*v)// bottom up
+int2d(Th3,299,301)(Bi*u*v)// Top down first outlet
+int2d(Th3,399,401)(Bi*u*v)// bottom up second outlet
+int2d(Th3,499,501)(Bi*u*v)// Top down second outlet
+int2d(Th3,599,601)(Bi*u*v)// bottom up third outlet
+int2d(Th3,699,700,701)(Bi*u*v)// Top down third outlet
+int3d(Th3)(-fh*v)
+int2d(Th3,100)(-qh*v)
;
//
// variational forms free from parameters. To be computed once
//
varf a1(u,v)= int3d(Th3,0,2,4,6,8,10)(Grad3(u)' *Grad3(v));// goes with factor k0
varf a2(u,v)= int3d(Th3,1,5,9)(Grad3(u)' *Grad3(v));// goes with factor kleft
varf a3(u,v)= int3d(Th3,3,7,11)(Grad3(u)' *Grad3(v));// goes with factor kright
varf a4(u,v)= int3d(Th3)( u*v ) ;// goes with factor aa
varf a5(u,v)=int2d(Th3,10,12,14,16,18,110,11,15,19,13,17,111
,30,32,34,36,38,310,31,35,39,33,37,311,
40,44,48,41,45,49,20,24,28,23,27,211,199,201
,299,301,399,401,499,501,599,601,699,700,701)
(u*v);// goes with factor Bi
//
// Right hand side also free from parameter. To be computed once
//
varf b1(unused,v)= int3d(Th3)( fh*v )+int2d(Th3,100)( qh*v );
// Matrices free from parameters out of the variational forms
//
matrix ma1= a1(Vh,Vh);
matrix ma2= a2(Vh,Vh);
matrix ma3= a3(Vh,Vh);
matrix ma4= a4(Vh,Vh);
matrix ma5= a5(Vh,Vh);
//
// Right hand side vector for problem
// Parameter free
//
// Computed once for all
//
Vh RHSunico;
RHSunico[]=b1(0,Vh);
//
// Matrices to gather parameter effects
//
//m1=k0*ma1;
//m2=kleft*ma2;
//m3=kright*ma3;
//m4=aa*ma4;
//m5=Bi*ma5;
//
matrix m1,m2,m3,m4,m5;
//
// Global matrix changes with parameter
//
matrix Aglobal;
//
//Aglobal=m1;
//Aglobal=Aglobal+m2;
//Aglobal=Aglobal+m3;
//Aglobal=Aglobal+m4;
//Aglobal=Aglobal+m5;
//
// Riesz Representation of RHS L(v)
// Parameter independent
//
include "RieszRepresentation_RHS.idp";
//
// SET A COLOR SCALE FOR DRAWING PLOTS (gracias Gladys)
//
real[int] colorhsv=[
4./6.,1, 0.5, //dark blue
4./6., 1, 1, //blue
0.5,1,1, //cyan
0.44, 1, 1, //light green
0.33,1,1, //green
0.2, 1, 1, //light yellow
0.15,1,1, //yellow
0.12,1,1, //orange
0.09, 1, 1, //light red
0.027, 1, 1 //red
];
//
// MAXIMUM SIZE OF REDUCED BASIS IS SAMPLES=20
//
int samples=20;
//
// FOR THE MASS MATRICES WITH REDUCED BASIS
//
// Ap(usample[i],usample[j]), p=0,1,2,3,4
//
// Remark: On settingANDproblem.idp I used notation a1,a2,a3,a4,a5
//
real[int,int] A0(samples,samples); //
real[int,int] A1(samples,samples);//
real[int,int] A2(samples,samples);//
real[int,int] A3(samples,samples);//
real[int,int] A4(samples,samples);//
//
// FOR THE RHS WITH THE REDUCED BASIS
//
real[int] Rrhs(samples); //VECTOR L(usample[i]) (0...9)//
//
real[int] cc(samples); // OUTPUTS
Vh[int] usample(samples); // REDUCE BASIS FUNCTIONS
real[int] mui(samples); // PARAMETERS
int[int] jbest(samples); // INDEX TO EACH PARAMETER
jbest=0;