Added matlab scripts for automatic femlab problem generation

git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@1838 f25d12d1-65f4-0310-ae8a-bbce733d8d8e

pde/DomainExternalSolverPackage/MFiles/femlabdomainfix.m

@@ -0,0 +1,53 @@
+function [neweb,poly]=femlabdomainfix(v,e);
+
+% v=[4.5 2 1;2 2 1;2 1 1;-0.5 1 1;-1.5 2 1;-4.5 2 1;-4.5 -2 2;-1 -2 3;-1 -1.5 3;-1.5 -1 3;-1.5 0 3;1.5 0 3;1.5 -1 3;1 -1.5 3;1 -2 3;4.5 -2 3];
+
+% e=[1 2 1;2 3 1;3 4 1;4 5 1;5 6 1;6 7 2;7 8 3;8 9 3;9 10 3;10 11 3;11 12 3;12 13 3;13 14 3;14 15 3;15 16 3;16 1 4];
+
+
+vx = v(:,1);
+vy = v(:,2);
+vb = v(:,3);
+
+es = e(:,1);
+ee = e(:,2);
+eb = e(:,3);
+
+poly=line2(vx,vy);
+
+newvertices=flgeomvtx(poly)';
+newedges=flgeomse(poly)';
+
+newvx=newvertices(:,1);
+newvy=newvertices(:,2);
+
+newes=newedges(:,1);
+newee=newedges(:,2);
+
+n=length(newes);
+neweb=-ones(n,1);
+
+myeps = eps*100;
+
+for i=1:n
+  news = newes(i);
+  newe = newee(i);
+  newsx = newvx(news);
+  newsy = newvy(news);
+  newex = newvx(newe);
+  newey = newvy(newe);
+  for j=1:n
+    s2 = es(j);
+    e2 = ee(j);
+    sx = vx(s2);
+    sy = vy(s2);
+    ex = vx(e2);
+    ey = vy(e2);
+    if ((abs(newsx-sx) < myeps) & (abs(newsy-sy) < myeps) & ... 
+        (abs(newex-ex) < myeps) & (abs(newey-ey) < myeps)) | ...
+       ((abs(newsx-ex) < myeps) & (abs(newsy-ey) < myeps) & ... 
+        (abs(newex-sx) < myeps) & (abs(newey-sy) < myeps))
+      neweb(i) = eb(j);
+    end
+  end
+end

pde/DomainExternalSolverPackage/MFiles/genfemlabheatproblem.m

@@ -0,0 +1,69 @@
+function fem=genfemlabproblem(filename,bcfilename,f,c,da,t2,N)
+clear fem;
+[v,e,vh]=readgeom(filename);
+[neweb,poly]=femlabdomainfix(v,e);
+
+[bctype,bcgval,bcqval,bcind]=readbc(bcfilename);
+%bcgroups=cell(size(bcind',1),1);
+%for i=1:size(neweb,1),
+%    bcgroups{neweb(i)}=[bcgroups{neweb(i)} i];
+%end
+fem.const={};
+for i=1:size(bcind',1),
+    if bctype(i) == 1 % dirichlet
+        gval{i}={{0}};
+        qval{i}={{0}};
+        hval{i}={{1}};
+        str=sprintf('bcrval_%d',i);
+        rval{i}={{str}};
+        fem.const=cat(2,fem.const, {str bcgval(i)});
+    elseif bctype(i) == 2 % neumann
+        gval{i}={{bcgval(i)}};
+        qval{i}={{0}}; % should be zero also in bcqval
+        hval{i}={{0}};
+        rval{i}={{0}};        
+    elseif bctype(i) == 3 % robin
+        str=sprintf('bcgval_%d',i);
+        gval{i}={{str}};
+        fem.const=cat(2,fem.const, {str bcgval(i)});
+        str=sprintf('bcqval_%d',i);
+        qval{i}={{str}};
+        fem.const=cat(2,fem.const, {str bcqval(i)});
+        hval{i}={{0}};
+        rval{i}={{0}};        
+    elseif bctype(i) == 4 % timedepdirichlet
+        warning('timedepdirichlet not implemented yet, using constant g value');
+        gval{i}={{0}};
+        qval{i}={{0}};
+        hval{i}={{1}};
+        %rstr=sprintf('''%f * t''',bcgval(i));
+        str=sprintf('timedepgval_%d',i);
+        rval{i}={{str}};        
+        fem.const=cat(2,fem.const, {str bcgval(i)});
+    end
+end
+%fem.dim='u';
+%fem.form='coefficient';
+fem.equ.f=f;
+fem.equ.init=0;
+fem.equ.c=c;
+fem.geom=poly; % see femlabdomainfix.m
+fem.bnd.ind=neweb; % see femlabdomainfix.m
+%fem.bnd.ind=bcgroups; % see femlabdomainfix.m
+fem.bnd.g=gval;
+fem.bnd.q=qval;
+fem.bnd.h=hval;
+fem.bnd.r=rval;
+fem.shape=2;
+hmax=max(vh); % Get the biggest h
+fem.equ.da=da;
+fem.mesh=meshinit(fem,'hmax',hmax,'Report','off');
+%fem.xmesh=meshextend(fem,'report','off');
+if da==0
+%    fem.sol=femlin(fem,'Report','off');
+else
+%    fem.init = asseminit(fem);
+%    fem.sol=femtime(fem,'tlist',linspace(0,t2,N));
+%    fem.sol=femtime(fem,'tlist',0:t2/N:t2);
+end
+

pde/DomainExternalSolverPackage/MFiles/genfemlabproblem.m

@@ -0,0 +1,62 @@
+function fem=genfemlabproblem(filename,bcfilename,f,c,da,t2,N)
+clear fem;
+[v,e,vh]=readgeom(filename);
+[neweb,poly]=femlabdomainfix(v,e);
+
+[bctype,bcgval,bcqval,bcind]=readbc(bcfilename);
+%bcgroups=cell(size(bcind',1),1);
+%for i=1:size(neweb,1),
+%    bcgroups{neweb(i)}=[bcgroups{neweb(i)} i];
+%end
+fem.const={};
+for i=1:size(bcind',1),
+    if bctype(i) == 1 % dirichlet
+        gval{i}={{0}};
+        qval{i}={{0}};
+        hval{i}={{1}};
+        rval{i}={{bcgval(i)}};
+    elseif bctype(i) == 2 % neumann
+        gval{i}={{bcgval(i)}};
+        qval{i}={{0}}; % should be zero also in bcqval
+        hval{i}={{0}};
+        rval{i}={{0}};        
+    elseif bctype(i) == 3 % robin
+        gval{i}={{bcgval(i)}};
+        qval{i}={{bcqval(i)}};
+        hval{i}={{0}};
+        rval{i}={{0}};        
+    elseif bctype(i) == 4 % timedepdirichlet
+        warning('timedepdirichlet not implemented yet, using constant g value');
+        gval{i}={{0}};
+        qval{i}={{0}};
+        hval{i}={{1}};
+        %rstr=sprintf('''%f * t''',bcgval(i));
+        rval{i}={{'timedepgval'}};        
+        fem.const=cat(2,fem.const, {'timedepgval' bcgval(i)});
+    end
+end
+%fem.dim='u';
+%fem.form='coefficient';
+fem.equ.f=f;
+fem.equ.init=0;
+fem.equ.c=c;
+fem.geom=poly; % see femlabdomainfix.m
+fem.bnd.ind=neweb; % see femlabdomainfix.m
+%fem.bnd.ind=bcgroups; % see femlabdomainfix.m
+fem.bnd.g=gval;
+fem.bnd.q=qval;
+fem.bnd.h=hval;
+fem.bnd.r=rval;
+fem.shape=2;
+hmax=max(vh); % Get the biggest h
+fem.equ.da=da;
+fem.mesh=meshinit(fem,'hmax',hmax,'Report','off');
+%fem.xmesh=meshextend(fem,'report','off');
+if da==0
+%    fem.sol=femlin(fem,'Report','off');
+else
+%    fem.init = asseminit(fem);
+%    fem.sol=femtime(fem,'tlist',linspace(0,t2,N));
+%    fem.sol=femtime(fem,'tlist',0:t2/N:t2);
+end
+

pde/DomainExternalSolverPackage/MFiles/genfemlabproblem2.m

@@ -0,0 +1,60 @@
+function fem=genfemlabproblem2(filename,bcfilename,f,c,da,t2,N)
+clear fem;
+[v,e,vh]=readgeom(filename);
+[neweb,poly]=femlabdomainfix(v,e);
+
+[bctype,bcgval,bcqval,bcind]=readbc(bcfilename);
+bcgroups=cell(size(bcind',1));
+%for i=1:size(eb',1),
+%    bcgroups{eb(i)}=[bcgroups{eb(i)} i];
+%end
+for i=1:size(bcind',1),
+    if bctype(i) == 1 % dirichlet
+        gval{i}=0;
+        qval{i}=0;
+        hval{i}=1;
+        rval{i}=bcgval(i);
+    elseif bctype(i) == 2 % neumann
+        gval{i}=bcgval(i);
+        qval{i}=0; % should be zero also in bcqval
+        hval{i}=0;
+        rval{i}=0;        
+    elseif bctype(i) == 3 % robin
+        gval{i}=bcgval(i);
+        qval{i}=bcqval(i);
+        hval{i}=0;
+        rval{i}=0;        
+    elseif bctype(i) == 4 % timedepdirichlet
+        warning('timedepdirichlet not implemented yet, using constant g value');
+        gval{i}=0;
+        qval{i}=0;
+        hval{i}=1;
+%        rstr=sprintf('''%f * t''',bcgval(i));
+        rval{i}=bcgval(i);        
+    end
+end
+%fem.dim='u';
+%fem.form='coefficient';
+fem.appl{1}.mode.class='FlPDEC';
+fem.appl.equ.f=f;
+fem.equ.init=0;
+fem.equ.c=c;
+fem.geom=poly; % see femlabdomainfix.m
+fem.bnd.ind=neweb; % see femlabdomainfix.m
+fem.bnd.g=gval;
+fem.bnd.q=qval;
+fem.bnd.h=hval;
+fem.bnd.r=rval;
+fem.shape=2;
+hmax=max(vh); % Get the biggest h
+fem.equ.da=da;
+fem.mesh=meshinit(fem,'hmax',hmax,'Report','off');
+%fem.xmesh=meshextend(fem,'report','off');
+if da==0
+%    fem.sol=femlin(fem,'Report','off');
+else
+%    fem.init = asseminit(fem);
+%    fem.sol=femtime(fem,'tlist',linspace(0,t2,N));
+%    fem.sol=femtime(fem,'tlist',0:t2/N:t2);
+end
+

pde/DomainExternalSolverPackage/MFiles/readbc.m

@@ -0,0 +1,18 @@
+function [bctype,bcgval,bcqval,bcind]=readbc(filename)
+fid=fopen(filename);
+n=fscanf(fid,' { %d , %d ',2);
+nbc=n(1);
+bcdim=n(2);
+if bcdim ~= 4
+    error('bcdim must be 4');
+end
+fscanf(fid,' , { ',1);
+for i=1:nbc,
+    val=fscanf(fid,' { %f , %f , %f , %f }', 4);
+    bctype(i)=val(1);
+    bcgval(i)=val(2);
+    bcqval(i)=val(3);
+    bcind(i)=val(4);
+    fscanf(fid,' , ',1);
+end
+fclose(fid);

pde/DomainExternalSolverPackage/MFiles/readgeom.m

@@ -0,0 +1,25 @@
+function [v,e,vh]=readgeom(filename)
+fid=fopen(filename);
+findToken(fid,'Vertices');
+nv=fscanf(fid,'%d',1);
+for i = 1:nv,
+    v(:,i)=fscanf(fid,'%f %f %d',3);
+end;
+findToken(fid,'Edges');
+ne=fscanf(fid,'%d',1);
+for i = 1:nv,
+    e(:,i)=fscanf(fid,'%d %d %d',3);
+end;
+findToken(fid,'hVertices');
+for i = 1:nv,
+    vh(i)=fscanf(fid,'%f',1);
+end;
+v=v';
+e=e';
+
+
+function findToken(fid, token)
+t=fscanf(fid, '%s', 1);
+while strcmp(token,t) == 0 || feof(fid)
+    t=fscanf(fid, '%s', 1);
+end;